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ABSTRACT 

Electromagnetic  scattering  from  trees  and  vegetation  is  of  prime  importance  in 
radar  and  remote  sensing.  The  actual  problem  of  scattering  from  trees  is  rather 
complicated  and  involves  three  dimensional  scattering  from  lossy,  electrically  large, 
and  randomly  oriented  objects. 

In  this  thesis,  the  radar  cross  section  of  a  planar  fractal  tree  is  considered. 
Although  a  planar  tree  is  far  from  being  real,  scattering  from  it  sheds  light  on  tlie 
scattering  phenomenon  from  an  actual  tree.  The  planar  tree  is  generated  using 
fractal  geometry  and  its  branches  are  considered  perfectly  conducting.  The  tree  is 
illuminated  by  a  plane  wave  and  the  problem  is  solved  using  the  moment  method. 
Data  is  presented  for  the  radar  cross  section  for  different  branching  angles  of  the 
tree  and  at  different  frequencies. 
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I.  INTRODUCTION 

A.   NEED  FOR  THE  STUDY 

The  trees  existing  in  the  natural  world  are  fractal  anisotropic.  They  are  made 
up  of  long,  intersecting  and  lossy  objects.  The  geometry  of  these  objects  is  not  easy 
to  set  up  as  they  are  randomly  oriented  in  a  three  dimensional  space.  The  use  of 
fractals  facilitates  the  modeling  of  semi— randomly  distributed  structures. 
Mandelbrodt  [Ref.  1]  has  shown  that  a  number  of  naturally  occurring  phenomenon 
such  as  coastlines,  clouds,  trees,  etc.  are  fractal  in  nature.  For  instance,  when  a 
branch  is  divided  into  two  (or  more),  the  ratio  of  the  length  of  the  subbranches  to 
the  main  branch  length  remains  constant.  Furthermore,  the  branching  angle  also 
remains  the  same. 

In  this  thesis,  scattering  from  planar  fractal  trees  is  considered.  The  radar 
cross  section  of  a  real  fractal  tree  is  a  complicated  scattering  problem.  The  analysis 
of  this  problem  in  a  two  dimensional  space  gives  an  approximate  idea  of  the  radar 
cross  section  of  a  real  tree. 

The  radar  cross  section  of  an  object  is  a  quantitative  measure  of  the  ratio  of 
the  power  density  that  is  received  and  scattered  by  the  object  to  the  power  density 
of  the  electromagnetic  wave  that  illuminates  that  object.  The  radar  cross  section  is 
independent  of  the  range  of  the  object  for  the  far-field  situation.  The  theoretical 
definition  of  the  radar  cross  section  "cr"  is  given  by  the  formula: 


where  e'  is  the  incident  electric  field  vector,  E  is  the  scattered  electric  field  vector, 
and  R  is  the  distance  between  the  scattered  object  and  the  point  of  observation.  The 
radar  cross  section  has  dimensions  of  area.  Usually,  it  is  expressed  in  square 
wavelengths. 

B.      STATEMENT  OF  THE  PROBLEM 

This  thesis  investigates  the  radar  cross  section  of  planar  fractal  trees.  These 
trees  are  composed  of  planar  thin  strip  dipoles  of  arbitrarily  dimensions  and 
orientations  in  a  plane.  These  structures  are  excited  by  plane  waves  of  various 
frequencies. 

The  first  step  in  solving  the  problem  is  to  calculate  the  induced  current 
distribution  on  each  strip.  The  calculation  of  the  current  distribution  is  based  on  the 
theory  of  the  moment  method  and  requires  a  knowledge  of  the  impedance  between 
any  two  of  these  strips  as  well  as  the  voltage  on  each  planar  strip  due  to  the 
incident  electric  field.  The  basic  concepts  and  the  calculation  of  the  current 
distribution  are  described  in  Chapter  2. 

In  Chapter  3,  the  development  of  a  FORTRAN  program,  is  presented.  The 
evaluation  of  the  radar  cross  section  requires  the  knowledge  of  the  scattered  electric 
field  due  to  the  induced  currents  on  the  planar  strips.  The  program  computes  the 
scattered  electric  field  and  then  the  radar  cross  section  of  that  structure.  The  details 
of  these  calculations  are  also  presented  in  this  Chapter. 

The  computer  models  that  are  used  by  the  developed  program  are  presented  in 
Chapter  4.  Their  generation  is  based  on  the  fractal  geometry.  An  existing  and 
modified  program  is  used  to  generate  the  geometry  of  the  planar  fractal  trees  in 
order  to  be  used  as  input  in  the  developed  program. 


The  numerical  results  of  the  radar  cross  section  of  a  single  planar  dipole  and  a 
a  number  of  planar  fractal  trees  are  presented  in  Chapter  5.  The  scattering  from  a 
single  dipole  is  compared  with  standard  results  for  a  similar  case.  The  limitations  of 
the  developed  program  are  also  presented  in  this  Chapter. 

In  Chapter  5,  the  conclusions  of  the  radar  cross  section  results  and 
recommendations  are  presented.  The  two  programs  that  are  used  for  this 
investigation  are  listed  in  the  Appendices. 


II.     METHOD  OF  MOMENTS  THEORY 

In  this  section  of  the  study,  the  basic  concepts  of  the  moment  method  theory 
are  presented.  This  theory  is  used  in  the  development  of  the  RCS  program  to  find 
the  current  distribution  on  a  planar  strip  due  to  an  incident  plane  wave. 

For  a  given  structure  consisting  of  planar  dipoles  the  impedance  between  any 
two  of  them  is  calculated  from  the  knowledge  of  the  geometry  and  the  wavelength  of 
the  incident  plane  wave.  The  voltage  on  each  dipole  is  calculated  from  the 
knowledge  of  the  characteristics  of  the  incident  electric  field.  The  induced  current 
distribution  on  each  dipole  of  the  structure  is  determined  from  the  calculated 
impedance  and  voltage  using  the  method  of  moments  theory. 

A.       GENERAL  THEORY 

The  method  of  moments  is  a  numerical  procedure  for  solving  integral 
equations  of  the  form: 

rb 

f(x')K(x,x')dx'  =  g(x),     a<x<b  (eqn  2.1) 

-^a 

where  f(x')  is  an  unknown  function,  K(x,x')  is  a  known  Kernel  or  Green's  function, 
and  g(x)  is  a  given  function.  This  procedure  reduces  the  integral  equation  (eqn  2.1) 
to  a  system  of  simultaneous  linear  algebraic  equations  in  terms  of  some  unknown 
coefficients.  This  method  requires  that  the  function  f(x')  be  approximated  by  a 
series  of  N  expansion  functions  or  "  basis  functions  ",  such  that 


f(x')^    !:anfn(x'),  n=l,2, ,N  (eqn  : 

n  =  l 


where  the  domain  of  fn(x')  is  the  same  as  that  of  f(x')  and  an's  are  the  complex 
unknown  expansion  coefficients. 

There  are  two  types  of  basis  functions.  The  subdomain  functions,  which  are 
nonzero  over  a  part  of  the  domain  of  the  unknown  function  f(x'),  and  entire  domain 
functions  being  nonzero  over  the  entire  domain  of  f(x').  In  antennas,  some 
commonly  employed  subdomain  basis  functions  are  the  piecewise  sinusoid  functions, 
the  unit  height— pulse  functions  and  the  piecewise  triangular  functions. 

The  subdomain  procedure  requires  subdivision  of  the  structure  into  N 
nonoverlapping  segments.  Figure  2.1  shows  a  segmented  line  where  the  segments  are 
assumed  to  be  coUinear  and  of  equal  length,  although  this  condition  is  not  necessary. 
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Figure  2.1   Segmented  Line.  [From  Ref.  2] 


Figure  2.2  shows  a  subdomain  unit  height— pulse  function  which  produces  a 
staircase  representation  of  the  unknown  function  f(x').  Figure  2.3  sliows  a 
subdomain  sinusoid  basis  function  and  the  representation  of  the  function  f(x'). 
Figure  2.4  shows  a  subdomain  triangular  basis  function  producing  a  smoother 
representation  of  the  function  f(x')  than  the  case  of  the  unit  height-pulse  basis 
function. 

The  use  of  entire— domain  basis  functions  does  not  require  any  segmentation  of 
the  structure.  One  of  the  most  most  commonly  used  basis  functions  of  this  kind  is 
the  sinusoidal  basis  functions. 


Figure  2.2   Unit  Height— Pulse  Basis  Function  and  a  Staircase 
Representation  of  f(x').  [From  Ref.  2] 


Figure  2.3   Sinusoid  Basis  Function  and  a  Representation 
off(x').  [FromRef.  3] 


Figure  2.4  Tiiangular  Basis  Function  and  a  Repiesentation 
off(x').  [FromRef.  2] 


The  substitution  of  the  function  f(x')  by  a  sequence  of  N  basis  functions  leads 
to  one  equation  of  N  unknowns  of  the  expansion  coefficients  an,  which  can  be  found 
by  using  N  linearly  independent  equations.  These  equations  are  set  up  by  the  use  of 
"testing  or  weighting"  functions  a;m(x),  m=l,2,....N.  At  this  point  the  definition  of 
the  inner  product  is  required.  The  inner  product  (f,  g>  between  any  two  functions  or 
vectors  f,  g  is  a  scalar  operation  defined  as 

<f,g>  =j|    f-gds  (eqn2.3) 

where  S  is  the  surface  of  the  structure  that  is  analyzed  [Ref  4].  The  inner  product  of 
the  selected  testing  functions  a;n)(x),  with  the  two  sides  of  the  original  integral 
equation  leads  to  the  equation: 


(^  a;4x),  j      f(x')K(x,xOdx'  ^  =  <(^g(x),  a;„,(x)  ^ 


(eqn  2.4) 


b    N 

m  =  1,2, ,N 

(eqn  2.5) 


<^a;m(x),  S  anfn(x')dxM  =  /g(x),  a;m(x)  y 


The  use  of  the  above  definition  for  the  inner  product,  yields: 


N      c    b 
San 
n=i    ^a 


r      U  r    U  r    D 

an      a;m(x)dx       fn(x')K(x,x')dx'  =        g(x)a)n,(x)dx,      m  =  1,2, N 

J  a.  J  a  -^  a 

(eqn  2.6) 
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The  following  substitutions 

pb 
\'m  =        g(x)^'4x)dx,         m=l,2, ,N  (eqn  2.7) 


and 


r  b  p  b 

Zmn  =        c^'m(x)dx        fn(x' )K(x,x' )dx' ,       m,n=l,2,....N 


(eqn  2.8) 


result  in  a  matrix  equation  of  the  form: 

[Zmn]-[anl  =  [Vj  (eqn  2.9) 

The  unknowns  [an]  can  be  obtained  through  a  matrix  inversion: 

lanl  =  [Zmnl''-[V„,]  (eqn  2.10) 

- 1 
where  [Zmn]      is  the  inverse  matrix.  The  column  vector  [Vm]  depends  upon  the  given 

function  g(x)  and  the  selected  testing  functions  a;n)(x).  The  matrix  [Zmn]  depends 

upon  the  known  kernel  K(x,x')  and  both  the  selected  basis  and  testing  functions. 

Once  the  expansion  coefficients  are  known,  the  function  f(x')  is  also  known. 

The  choice  of  basis  and  testing  or  weighting  functions  is  based  upon  experience 

and  the  rule  is  that  their  number  has  to  be  the  same.  The  procedure  of  using  the 

same  basis  and  testing  functions  is  called  the  "  Garlekin's  method  ". 

B.      APPLICATIONS  TO  EM  THEORY 

A  number  of  problems  in  electromagnetic  radiation  and  scattering  can  be 
solved  by  the  method  of  moments.  In  the  present  case,  the  simple  case  of  a  perfectly 
conducting  object  "  situated  in  free  space  "  is  considered. 

The  basic  problem  is  to  investigate  the  case  when  an  object  is  illuminated  by 
fields  of  known  impressed  electric  and  magnetic  currents  (j\  M^  In  the  absence  of 


the  object,  the  impressed  currents  radiate  the  assumed  known  incident  electric  and 


radiate  the  unknown  total  fields  (E  ,  H  ). 

The   integral   equation   is   obtained   by   the   surface   equivalence   principle, 
replacing  the  object  by  free  space  together  with  the  electric  surface  current  density 

J  =  nxH^  (eqn2.11) 

where  J  exists  on  the  entire  surface  S  of  the  object  and  n  is  the  unit  vector  normal 


E^  =  E^-E^  (eqn2.12) 

H^  =  H^-H^  (eqn2.13) 

The  boundary  conditions  for  this  case  enforce  the  total  tangential  electric  field 
on  the  surface  S  to  zero: 

n  X  (  E^  4-  E^  )  =  0  (eqn  2.14) 

This  is  an  integral  equation  for  J  since  the  sc.ittered  electric  field  E  can  be  written 
as  an  integral  over  S  of  the  dot  product  of  J  and  the  dyadic  free  space  Green's 
function.  [Ref.  5] 

As  the  geometry  of  the  object  is  known,  it  is  more  convenient  to  use  the  total 
current  I  instead  of  the  total  current  density  J.  So,  the  problem  is  to  find  the 
unknown  current  I  induced  by  the  incident  electric  field.  The  method  of  moments 
solve  these  kind  of  problems  by  the  following  procedure. 

The  first  step  is  to  expand  the  unknown  current  I  in  terms  of  some  basis  set: 

N 
I«    SInFn  (eqn  2.15) 

n  =  l 

where  the  In  are  the  sequence  of  N  unknown  complex  coefficients,  and  the  Fn  is  a 
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sequence  of  N  known  modes  or  basis  functions.  The  best  choice  of  Fq  for  a  given 
problem  could  be  quite  involved  and  is  discussed  in  [Ref.  4,  pp  308—310]. 

The  second  step  is  to  select  the  testing  or  weighting  functions  a;ni,  m=l,2,....N. 
These  can  be  identical  with  the  basis  functions  or  different.  The  inner  product  of  the 
sequence  of  N  weighting  functions  Um,  with  both  sides  of  the  integral  equation  gives 
a  N  X  N  system  of  simultaneous  linear  algebraic  equations  of  the  symbolic  form: 
[Z].[I]  =  [V]  (eqn2.16) 

where  I  is  the  current  column  vector  whose  N  components  give  the  values  of  In 
and  [Zj  is  the  N  x  N  impedance  matrix  given  by  the  equation 

Zmn  =  -       En-o^nids  (equ  2.17) 

where  S  is  the  surface  of  the  structure  being  analyzed.  The  impedance  matrix  [Z]  is 
always  symmetric,  and,  for  the  special  case  of  a  thin  dipole  instead  of  an  arbitrary 
object,  is  also  a  Toeplitz  matrix.  In  a  Toeplitz  matrix,  Zmn  depends  only  on  |m— n  | . 

Generally,  [Z]  is  dependent  only  on  the  geometry  and  material  composition  of 
the  scatterer,  but  not  on  the  incident  fields  [Ref  5].  The  right-hand  side  of  the  last 
equation  is  the  voltage  vector  whose  N  components  give  the  corresponding  mode 
voltage.  The  voltage  vector  depends  only  on  the  excitation,  i.e.,  the  incident  electric 
field.  The  dimensions  of  the  elements  of  [Z]  and  [V]  are  volt-amps  (VA),  while  the 
elements  of  [I]  are  dimensionless. 

The  solution  of  the  last  matrix  equation  is  the  column  vector  [I],  whose 
elements  represent  the  complex  coefficients  In-  As  the  total  current  I  was  expanded 
by  a  sequence  of  N  known  modes  or  basis  functions  and  the  complex  coefficients  In 
are  known,  the  total  current  and  so  the  total  current  density  J  is  also  known. 
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Although  the  choice  of  weighting  functions  is  free,  it  has  to  be  considered  that 

2 
the  matrix  equation  being  solved  requires  the  evaluation  of  N    terms  and  each  term 

requires   two   or   more   integrations.    When   these   integrations   are   to   be   done 

numerically,  the  computations  become  complicated.  There  is  a  way  to  reduce  this 

complexity  by  choosing  as  weighting  functions  the  Dirac  delta  functions.  This  is  the 

method  of  point -matching  in  which  delta  functions  are  enforced  only  at  discrete 

points  on  the  surface  S.  The  results  of  this  method  can  be  quite  accurate  especially 

when  the  discrete  points  are  selected  to  be  equally  spaced.  The  solution  satisfies  the 

electromagnetic  boundary  conditions  (e.g.,  vanishing  tangential  electric  fields  on  the 

surface  of  an  electric  conductor)  only  at  discrete  points  [Ref  4].  Between  these 

points  the  boundary  conditions  may  not  be  satisfied.  In  this  case  it  is  required  to 

define  the  deviation  as  a  residual  (e.g.,  residual=AE|tan  =  E  |tan+  E^itan  /  0  on 

the  surface  S  of  an  electric  conductor)  and  use  the  method  of  weighted  residuals  so 

that  the  boundary  conditic^ns  will  be  satisfied  in  an  average  sense  over  the  entire 

surface  S. 
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III.    ANALYSIS  AND  DEVELOPMENT  OF  RCS  PROGRAM 

In  this  chapter,  the  basic  steps  that  are  involved  in  the  calculation  of  the 
Radar  Cross  Section  (RCS)  of  a  planar  structure  composed  of  conducting  strips  are 
presented.  As  mentioned  earlier,  the  fractal  tree  is  modeled  as  consisting  of  planar 
thin  strips.  The  radar  cross  section  of  the  fractal  tree  is  calculated  using  the  method 
of  moments.  A  Fortran  program  is  developed  to  calculate  the  RCS  of  the  tree  for  a 
specified  geometry  and  at  a  given  frequency. 

The  moment  method  discretizes  an  integral  equation  to  a  matrix  equation  of 
the  form: 

[Z]-[I]  =  [V]  (eqn.3.1) 

where  [Z]  is  the  impedance  matrix  whose  elements  represent  the  mutual  impedance 
between  any  two  dipoles  of  a  given  structure  depending  upon  the  wavelength  and 
geometry  of  the  structure,  [V]  is  the  voltage  matrix  whose  elements  correspond  to 
the  voltage  on  each  dipole  due  to  excitation  (  incident  electric  field  ),  and  [I]  is  the 
unknown  matrix  whose  elements  represent  the  induced  current  on  each  dipole. 

A.      DEVELOPMENT  OF  RCS  PROGRAM 

The  numerical  results  for  RCS  are  obtained  from  a  Fortran  program  that 
calculates  backscattered  RCS  of  a  planar  structure  consisting  of  arbitrarily  oriented 
dipoles.  For  convenience,  the  structure  will  be  assumed  to  be  in  the  y-z  plane  of  an 
xyz  cartesian  coordinate  system. 

Figure  3.1  shows  the  structure  to  be  investigated.  A  large  number  of  planar 
dipoles    of    variable    lengths,    widths,    and    orientations    in    the    y— z    plane    is 
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illuminated  by  a  plane  wave  with  electric  field  linearly  polarized  and  characterized 

by  the  angles  (/?o  and  ^o-  It  is  required  to  calculate  the  backscattered  RCS  from  this 

structure. 

The  incident  electric  field  E^  is  given  by  the  formula 

p='^.^-K^^-^  +  ^ry  +  ^^'^)  (eqn3.2) 

where  t  is  the  propagation  vector  and 

kx^+  ky^+  kz^  =  ko^  =  u?HQ(o  (eqn  3.3) 

ko  =  -^^  (eqn  3.4) 

A 

kx  =  kosin^ocoscpo  (eqn  3.5) 

ky  =  kosin^osintpo  (eqn  3.6) 

kz  =  kocos^o  (eqn  3.7) 

e  =  Exx  +  Eyy  +  EzZ  (eqn  3.8) 

The  electric  field  E^  lies  on  the  plane  perpendicular  to  the  direction  of 

propagation  of  the  plane  wave.  Hence,  e-k  =  0  implies  that 

Exkx  +  Eyky  +  Ezkz  =  0  (eqn  3.9) 

The  substitution  of  Ey  and  Ez  by  the  variables  a  and  b  (independent  variables)  gives 

the  coordinate  Ex  as  a  dependent  variable 


g^^_(ak^^_+^_bkzl  (eqn  3.10) 

kx 


Equation  3.1  the  becomes 

fii  =  [  ay  +  bz-(  ^^v  +  ^^^  K]e-'^^  ^-"^  +  ^^^  +  ^^'  ^ 
^  kx  J 

(eqn  3.11) 
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es-c^ 


Figure  3.1   Geometry  of  RCS  Program. 
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The  RCS  program  that  has  been  developed  makes  the  following  assumptions: 

1.  The  dipoles  are  very  thin  planar  strips.  The  width  of  the  dipoles  is 
assumed  to  be  electrically  small  so  that  only  the  axial  current  is 
significant.  Further,  there  is  no  current  variation  along  the  width  of 
the  dipoles. 

2.  The  dipoles  are  not  intersecting. 

3.  The  dipoles  are  perfect  conductors. 

The  program  uses  the  theory  of  moment  method  and  the  way  that  it  solves  the 
problem  can  be  understood  if  a  single  planar  thin  dipole  in  the  y-z  plane  is 
considered.  Each  dipole  is  subdivided  into  equal  segments.  Overlapping  Piecewise 
Sinusoidal  (PWS)  modes  are  assumed  to  exist  on  the  dipole.  The  length  of  each 
PWS  mode  is  equal  to  the  length  of  two  segments.  Figure  3.2  shows  the  nith  PWS 
mode.  The  PWS  mode  has  a  length  2hn,,  and  a  width  2wn).  The  coordinates  of  its 
center  are  (ym,  Zm),  and  it  is  oriented  at  an  angle  ipm  measured  from  the  z  axis. 

The  incident  electric  field  induces  current  along  the  axis  (  of  that  mode.  In 
Figure  3.2,  the  induced  current  J  varies  sinusoidally  along  the  axis  of  the  mth  mode. 
A  new  coordinate  system  r/— (  is  introduced,  where  (  is  the  axis  of  the  dipole  and  t]  is 
the  axis  perpendicular  to  (.  The  ?/— (  coordinate  system  is  obtained  by  rotating  the 
y-z  system  about  the  x-axis  through  an  angle  90^—ipjn-  Figure  3.3  shows  the  details 
of  this  rotation. 

The  relation  between  the  coordinates  of  the  center  of  the  mode  along  the  axis 
of  the  two  orthogonal  systems  is  found  to  be: 

y  =  yn,  +  Csin(Vn.)  -  ^os(Vm)  (eqn  3.12) 

z  =  Zn,  +  Ccos(V'm)  +  ^in(Vn))  (eqn  3.13) 
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The  corresponding  vector  equation  is: 

y  =  Csin(t'ni)  -  r;cos(V'm) 

Z  =  CCOS(t;m)  +  mHtm) 


(eqn3.14) 
(eqn  3.15) 


Fi«^ure  3.2   Geometry  of  the  mth  PWS  Mode  of  the  Structure. 


Figure  3.3  Transformation  of  PWS  Mode's  Center  Coordii 
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1.  Voltage  Equations 

In  this  section,  expressions  for  the  elements  of  the  voltage  matrix  [V] 
due  to  the  incident  electric  field,  are  presented.  Each  element  of  this  matrix 
corresponds  to  a  PWS  mode  of  the  structure  that  is  investigated.  The  size  of  the 
voltage  matrix  is  equal  to  the  total  number  of  modes  of  the  structure. 

The  induced  current  density  on  the  mth  mode  of  the  structure,  due  to 
incident  electric  field,  is 


J^^sin(ko(h.  -  KD)  (eqn3.16) 

2wn,      sin(kohn,) 


where  2\Vn,  is  the  width  and  2hni  is  the  length  of  the  mth  mode.  For  each  mode,  the 
induced  current  variation  will  be  sinusoidal  along  its  axis,  being  maximum  at  the 
center  and  zero  at  the  end  points.  The  corresponding  voltage  Vm  is 

r    a;m    r    hm  ^l 

V     =  E      •    JdCd//  (eqn3.17) 

"^      J-L.'J-h.     |x=0 

From  equations  3.11  and  3.16,  it  is  seen  that 

E'-J  =  [a  sin(V^,)  +  b  cos(i/;,)]  e-'A^y^^-  +  ^^*"(^^)  +  ^^('-^  +Ccos(^.)] 

(eqn  3.18) 
where  rpja  is  the  angle  between  the  z  axis  (reference  for  angle  measurements)  and  the 
axis  of  the  mode  and  ym,  Zm  are  the  coordinates  of  the  center  of  the  mode  along  the 
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y  and  z  axis  respectively.  The  substitution  into  equation  3.17,  gives  the  following 
form  of  Vm: 

e-j(kyym   +    kzZm)  r    K 

""  "       sin(koh„)        ''  '"''*"'  ^  "  ~'*'*"^'  J-h„  ^"^'^^  si„(k„(h„-|  CI  ))dC 

(eqn  3.19) 
where  k^  =  kySin(V^tn)  +  kzCOs(Vm)-  A  closed  form  evaluation  of  the  integral  in  this 
equation  leads  to  final  form  of  Vm 

Vn,  =     2^'^^'°"'2     [  cos(kohn,)  -  cos(k^h„,)  ]  k^  ^  ±  ko 

k^    -  ko 

(eqn  3.20) 

Vm  =  Vom  hm  sin(kohm)  k^  =  ±  ko 


(eqn  3.21) 


where 


e-j(kyym   +    kzZm)    r 

Vom  = a  sin(0m)  +  b  cos(^m) 

sin(kohm)         ^ 


(eqn  3.22; 


2.  Radiation  Equations 

In  this  section,  expressions  for  the  far-^one  scattered  fields  due  to  the 
mth  PWS  mode  are  developed.  For  far-field  observations  the  electric  field  is  given 
in  spherical  coordinates  by  the  following  equations  [Ref  4]: 
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Ej.  s:  0  (eqn  3.23) 

E  .  ^  -  ^^ (  L,^  +  r;  N  .  )  (eqn  3.24) 

E^  ^  +  ^^ (  L^  +  t;  N^  )  (eqn  3.25) 

where 

l^=  U    [M^cos^os;)+Mycosfcin^-Mzsin^e"^j^o^'^^^^'"'ds' 

(eqn  3.26) 

(eqn  3.27) 
N^  z=  U    [J:,cos^coSyJ  +  Jycosfein^-  JzSin^]  e+J^or'cosV'm  ^g/ 

(eqn  3.28) 

(eqn  3.29) 
T]=  120-  (eqn  3.30) 

The  quantities  Jx,  Jy.  Jz,  are  the  components  of  the  electric  current 
density  Jg  that  are  induced  on  the  mth  mode  over  the  surface  S,  and  the  quantities 
Mx,  My,  Mz.  are  the  coordinates  of  the  magnetic  current  density  Ms  over  the  surface 
S.  As  the  structure  is  on  the  y— z  plane  and  Ms  is  zero,  the  quantities  J  ,  L^,  L  are 
zero.  Equations  3.24  and  3.2.5  become: 


E    ::  3i]ioe ijZx  (eqn3.31) 

"470-'^ 


E    ^  ±ikoLil!!^2N  (eqn  3.-32) 

^  47rr  ^ 
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As  the  current  density  Jm  is  along  the  (-axis,  all  quantities  within  the 
integrals  for  J  and  M  will  be  expressed  in  terms  of  the  r^(  system.  The  surface 
element  that  is  used  in  all  integrals  is  ds'  =  dydz  =  d^dr;  =  d(,  as  the  width  of  the 
mode  is  assumed  to  be  very  small  in  terms  of  its  length.  The  transformations  from 
the  y— z  system  to  the  (—r)  system  are  made  by  using  the  following  substitutions: 

r'cos(V'm)  =  ys'm9s'm(p  +  zcos^  (eqn  3.33) 

where 

y  =  ym  +  Csin(V'm)  (eqn  3.34) 

z  =  Zm  4-  Ccos(V'm)  (eqn  3.35) 

and 


Jm  =  Jyy  +  Jzz  =  C  [ sin(ko(h„,-|  CI ))] 

L    2wn,sin(kohn,)  ^ 


(eqn  3.36) 


where 


Jy  =  J^  sin(V'm)  (eqn  3.37~ 

Jz  =  J^  cos(^m)  (eqn  3.38) 

These  substitutions  and  algebraic  manipulations  give  the  quantities  N^_^^  and  N 

for  the  mth  mode  in  the  (—7/  system: 


N^-N, 


-2ko 


^Om  -  '^$0 2        2    [  cos(En,hn,)  -  cos(kohn,)  j         Eni  ^  ±ko 

Em  -  ko 

(eqn  3.39) 

(eqn  3.40) 


and 
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^^l  ^  ^y^O 2^^^   2     [  COS(En,h,n)  -  COS(kohn,)  j  En,  /  ±ko 

Em   -   ko 

(eqn  3.41) 

(eqn  3.42) 
where 


1^'     ^  sini'mcosfein;^  -  cosi^'mSin^    j    J[ko(ymSin^in(^+Zn,cos^)] 
'  ^0  sin(kohn>) 

(eqn  3.43) 
V      _  siiu)niC0Sj2_  T    J[ko(ymSin^in(^+Zn,cos^)] 
^^0      sin(kohn,) 

(eqn  3.44) 
Em  —  ko(sinVniSinfein(;)  +  cos^mCOS^) 

(eqn  3.45) 
In    this    thesis,   only   the   monostatic   radar   cross    section    will    be 
considered.  In  the  radiation  equations  the  angles  $  and  ^p  represent  the  orientation 
angles  of  the  scattered  electric  field  E  .  Their  relation  with  the  incident  angles  ^o 
and  ^0  for  the  monostatic  case  is: 

e=z-eo  (eqn  3.46) 

ip  =  TT  —  v^o  (eqn  3.47) 

If  M  denotes  the  total  number  of  PWS  modes  in  the  structure  under 
investigation,  then 

m=l 

N^=    S    N  (eqn  3.49) 

"^     m-1     '^ 
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The  final  expression  for  the  equations  3.24  and  3.25  is 


^d=^'     ^  ^9m  (eqn3.50) 


where 


47rr 


3.  RCS  Equations 

When  an  incident  plane  wave  with  electric  field  E   strikes  the  object 
and  E^  is  the  scattered  electric  field,  the  radar  cross  section  a  is  defined  as 


a^  lim47rr^     1^  1^  (eqn  3.53) 


|EM 


For  the  scattered  electric  field  E",  its  spherical  coordinates  E^,  and  E  have  been 
calculated  by  the  equations  3.50  and  3.51.  The  incident  electric  field  E^  is  known  by 
means  of  equation  3.11. 


•  2  2  2       I  |2 

|E^|     =  |a|     +  |b|     +    _§:kx^L_bk^  (eqn  3.54) 

'  kx         1 


|E^|    =|E^|    +|E^|     =|C|     [[|N^|    +|N^| 
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(eqn  3.55) 


The  combination  of  equations  3.11,  3.50-3.52,  3.54-3.55,  and  3.53  leads  to  the  final 
formula  for  RCS: 

^^  ^a|     +    |b|     +       ^v  +  ^^- 

(eqn  3.56) 

This  is  the  formula  that  the  program  uses  to  compute  the  RCS  of  a 
planar  structure  for  a  given  set  of  incident  angles  ^o,  and  v?o- 

B.      INPUT-OUTPUT  OF  RCS  PROGRAM 

The  RCS  program,  which  is  described  in  Appendix  A,  is  a  FORTRAN 
program  having  two  input  files.  The  first  input  file,  INPUTl,  contains  the  data  that 
characterize  the  incident  plane  wave  and  the  data  that  describe  the  geometry  of  the 
planar  structure  whose  broadside  RCS  is  measured.  The  second  input  file,  INPUT, 
contains  the  set  of  incident  angles  ^o,  Vo- 

The  program  reads  from  the  INPUTl  file  the  following  input  data: 

1.  frequency  of  the  incident  plane  wave  in  GHz, 

2.  parameters  a  and  b  that  characterize  the  polarization  of  the  incident 
electric  field. 

3.  number  of  dipoles  that  the  structure  consists  of, 

4.  half  length  of  each  dipole  in  cm, 

5.  half  width  of  each  dipole  in  cm, 

6.  coordinates  in  cm  of  the  center  of  each  dipole  along  the  two  axes  of 
the  orthogonal  system  that  the  structure  lies, 
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7.  orientation    angle    of    each    dipole,    measured    from    the    vertical    axis, 
positive  in  the  clockwise  direction,  and 

8.  number  of  segments  of  each  dipole. 

As  the  program  reads  these  input  data  it  generates  the  geometry  of  the  PWS 
modes,  calculates  the  impedance  elements  between  any  two  modes  of  the  structure, 
and  fills  the  impedance  matrix  [Z].  The  size  of  the  impedance  matrix  depends  on  the 
total  number  of  modes.  When  this  matrix  is  calculated  and  filled,  the  program 
inverts  it  to  [Z]~   and  stores  it  for  later  use. 

The  next  step  is  to  read  from  the  INPUT  file  the  sets  of  incident  angles  ^o,  '■Po 
and  for  each  one  it  calculates  the  voltage  on  each  mode,  filling  the  voltage  matrix 
[V].  This  matrix  is  a  column  vector  which  depends  upon  the  excitation  only.  Then  it 
multiplies  the  stored  inverted  impedance  matrix  [Z)~"  by  the  voltage  matrix.  This 
yields  the  current  vector  [I]: 

[I]  =  [Z-^-[V]  (eqn3.57) 

As  the  induced  currents  are  known,  the  program  uses  the  previously  described 
radiation  equations  and  calculates  the  RCS  of  the  structure  corresponding  to  the 

given  set  of  incident  angles  ^o,  V'o-  The  output  RCS  is  normalized  to  the  square  of 

2 
the  wavelength  A  ,  and  is  given  in  dB.  The  program  reads  the  next  set  of  incident 

angles  ^o,  '/'o  and  repeats  the  same  procedure  to  compute  the  RCS  of  the  new  set. 

Although  the  input  lengths  and  widths  are  in  cm,  the  program  considers  them 

normalized   to   the   wavelength.    For   accurate   results   at   least   4   segments   per 

wavelength  are  chosen.  The  selection  of  the  width  of  each  dipole  is  arbitrarily  taken 

as  L/W  =  33. 
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IV.     FRACTAL  TREE  GENERATION 

The  problem  of  measuring  the  radar  cross  section  of  a  real  natural  tree  is 
complicated  as  it  requires  an  investigation  in  three  dimensional  space  with  very 
long,  intersecting  elongated  objects  that  are  randomly  oriented. 

For  simplicity  the  radar  cross  section  (RCS)  problem  is  solved  in  a  two 
dimensional  space.  RCS  is  calculated  from  a  planar  fractal  tree  whose  branches  are 
considered  as  thin  planar  dipoles  of  different  lengths  and  widths.  In  an  actual  tree, 
the  branches  are  lossy  and  in  general  anisotropic.  However,  in  the  present  model, 
the  tree  is  comprised  of  perfectly  conducting  branches.  The  geometry  of  this  tree  is 
based  upon  the  fractal  geometry. 

A.       DEFINITIONS 

Fractal  is  a  mathematical  set  or  object  whose  form  is  extremely  irregular  or 
fragmented  at  all  scales  [Ref.  6].  The  requirement  to  describe  the  shape  of  many 
objects  that  appear  in  the  natural  world,  such  as  trees,  mountains,  coastlines,  etc., 
led  to  the  generation  of  fractal  geometry.  As  Euclidean  geometry  cannot  give 
mathematical  expressions  to  describe  fractal  objects,  fractal  geometry  is  used  to 
describe  mathematically  many  natural  patterns. 
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The  basic  characteristics  of  fractal  objects  are  [Ref.  7]: 

1.  A  large  degree  of  heterogeneity. 

2.  A  self-similar  structure  over  many  size  scales.  Self-similarity  refers 
to  the  general  preservation  of  form  or  characteristic  regardless  of  the 
scale  of  observation. 

3.  The  lack  of  a  well-defined  (characteristic)  scale. 

The  geometric  characteristics  of  fractal  objects  are  useful  for  describing  phenomena 
of  nature,  such  as  scattering  of  objects  Hke  landscapes  or  surface  cracks.  Although 
these  objects  are  irregular  and  randomly  oriented  in  nature,  they  show  structural 
similarities  on  several  different  discrete  size  scales  [Ref.  7]. 

One  measure  of  structural  complexity  is  the  fractal  dimension  Df.  There  are 
several  definitions  of  Df  depending  upon  the  particular  application.  For  the  case  of 
fractal  trees  the  fractal  dimension  Df  is  the  measure  of  the  space  that  a 

self— similar  structure  fil.s,  and  it  varies  with  the  branching  levels  that  the  structure 
consists  of. 

The  most  useful  terms  from  fractal  geometry  that  are  used  to  describe  a  planar 
fractal  tree  are  the  following: 

1.  The  number  of  branch  segments  N  formed  from  each  preceding  branch 
segment. 

2.  The  constant  similarity  ratio  "r"  that  relates  the  fractional  reduction 
in  segment  length  for  each  segment  to  previous  level,  this  factor  is  less 
than  1.0. 

In  this  case,  the  fractal  dimension  Df  is  given  by  the  formula  [Ref.  7]: 

,  _  j^j,  _    log(total  branch  length)      _       log(rN) 
log(average  branch  length)  log(l/r) 
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This  equation  is  accurate  for  symmetric  structures  or  asymmetric  structures  having 
a  normal  distribution  function  [Ref.  7].  More  details  about  fractal  geometry  are 
discussed  in  [Ref.  1]  and  [Ref.  7]. 

B.      FRACTAL  TREE  GENERATION 

In  the  planar  fractal  structures  that  are  used  for  calculating  the  radar  cross 
section,  the  number  of  branch  segments  N  formed  from  each  preceding  branch 
segment  is  2.  The  principal  properties  of  these  trees  are  that  the  branches  are  not 
overlapping,  and  the  angle  $  between  any  two  branches  is  the  same.  The 
nonoverlapping  between  the  branches  of  each  fractal  model  is  achieved  by  choosing 
the  branch  angle  6  for  a  given  reduction  factor  by  the  following  empirical  formula 
[Ref.  7] 


Minimum  6  =  32.34  x  r^""^^  (Radians) 


This  angle  is  given  in  radians  and  r  is  the  desired  reduction  factor  which  is  a 
constant  for  the  structure.  As  the  reduction  factor  r  increases  the  tree  fills  more 
space. 

Five  types  of  planar  fractal  trees  are  considered  in  this  thesis.  Each  type 
corresponds  to  a  set  of  values  for  reduction  factor  r  and  minimum  branch  angle  0. 
Table  4.1  shows  these  sets  of  r  and  minimum  9,  where  9  is  given  in  degrees. 

Each  model  is  generated  by  selecting  the  desirable  reduction  factor  r  (r  <  1.0), 
the  corresponding  branch  angle  $  given  by  equation  4.1,  the  initial  length  from 
which  the  reduction  will  start,  the  value  of  N  (which  in  the  investigated  case  is  2), 


and  the  number  of  the  desired  levels  for  branch  generation.  Each  level  contains  2^ 
points  corresponding  to  the  end  points  of  the  reduced  length  branches. 


TABLE  4.1 
NUMERICAL  VALUES  OF  r  AND  MINIMUM  ANGLE  6. 

REDUCTION  FACTOR  fr)  ANGLE  d 
0.53  29.940 

0.55  46.430 
0.60  95.890 

0.66  145.350 

0.71  180.000 


Figures  4.1^.5  show  the  planar  fractal  trees  that  correspond  to  each  set  of  r 
and  6  of  Table  4.1  respectively.  As  the  branch  angle  6  increases,  the  spreading  of  the 
branches  becomes  larger.  In  all  trees  the  physical  length  of  the  initial  branch  is 
chosen  as  two  centimeters.  These  models  were  generated  by  a  FORTRAN  program 
which  has  the  following  input  data: 

1.  The  initial  length.  This  the  length  of  the  first  dipole  whose  length  will 
be  reduced  by  the  constant  reduction  factor  (r). 

2.  The  value  of  N,  being  2  in  all  cases. 

3.  The  initial  point  that  the  reduction  starts.  This  point  is  the  end  point 
of  the  initial  wire. 
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4.  The  number  of  desirable  levels  for  branch  generation. 

5.  The  branch  angle  6. 

The  program  generates  two  output  files.  The  first  is  the  input  file  for  the  RCS 
program  containing  all  the  required  data  that  were  mentioned  in  Chapter  3.  The 
second  output  file  contains  the  coordinates  of  the  start  and  end  points  of  each 
generated  branch. 

For  the  models  characterized  by  reduction  factors  0.53,  0.55,  and  0.60,  the 
program  limits  the  structure  size  when  the  length  of  a  branch  becomes  smaller  than 
0.1  of  the  wavelength  of  the  incident  plane  wave.  For  the  models  characterized  by 
reduction  factor  0.66  and  0.71,  this  limit  is  taken  as  0.15.  The  reason  is  that  as  the 
reduction  factor  increases  the  number  of  branches  in  a  particular  branch  level 
increases.  The  size  of  the  tree  is  truncated  so  that  the  number  of  unknowns  is 
manageable. 

The  geometry  of  all  planar  fractal  trees  is  in  the  same  coordinate  system  that 
the  RCS  program  uses.  All  models  have  been  generated  in  the  y-z  plane  of  a 
cartesian  coordinate  system. 
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Figure  4.1   Fractal  Tree  for  r  =  0.53  and  6  =  29.^ 
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Figure  4.2  Fractal  Tree  for  r  =  0.55  and  6  =  46.430. 
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Figure  4.3   Fractal  Tree  for  r  =  0.60  and  6  =  95.89°. 
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Figme  4.4    Fractal  Tree  for  r  =  0.6G  and  0  =  145.35^. 
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Figure  4.5  Fractal  Tree  for  r  =  0.71  and  ^  =  180°. 
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V.    NUMERICAL  RESULTS 

In  this  chapter,  computed  results  for  the  radar  cross  section  a  of  planar 
structures  are  presented.  In  order  to  check  the  program,  the  radar  cross  section 
results  of  a  single  planar  dipole  are  compared  with  the  results  given  in  [Ref.  8]. 
Details  of  this  comparison  are  presented  in  the  following  section. 

In  all  models  where  numerical  results  for  the  radar  cross  section  are  presented, 
the  following  factors  have  been  considered: 

1.  E-plane  is  the  plane  corresponding  to  <^o  =  QO,  and  ^o  varying  from  QO  to 
1800. 

2.  H— plane  is  the  plane  corresponding  to  6o  =  90^,  and  ipo  varying  from  —90^ 
to  +900. 

3.  The  resulting  RCS  a  is  normalized  to  the  square  of  the  wavelength  A  of  the 
incident  plane  wave  (<7/A2). 

4.  The  number  of  segments,  that  each  dipole  is  subdivided,  is  taken  as  four 
per  wavelength. 

5.  The  physical  dimensions  and  orientations  of  the  dipoles  used  to  construct  a 
planar  structure  are  the  same  when  this  structure  is  investigated  at  different 
frequencies,  and  only  the  segmentation  is  different  depending  upon  the  wavelength  A 
of  the  incident  plane  wave. 

6.  Each  PWS  mode  has  length  equal  to  the  length  of  two  segments. 

7.  The  incident  electric  field  is  linearly  polarized  and  the  parameters  a  and  b 
are  taken  as  0  and  1  respectively  for  this  investigation.  This  corresponds  to  having 
onlv  a  horizontal  magnetic  field. 
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A.      SCATTERING  FROM  A  SINGLE  DIPOLE 

In  Figure  5.1  a  centered  loaded  vertical  planar  thin  dipole  of  length  L  and 
width  2w  is  oriented  along  the  z  axis.  This  dipole  is  excited  by  a  plane  wave  of 
frequency  30  GHz  traveling  in  the  x-y  plane  (^o  =  90^,  (po  =  QO).  The  same 
situation  is  described  in  [Ref.  8,  pp.  510—515]  for  a  cylindrical  dipole  of  radius  a  and 
length  L  such  that  L/2a=74.2.  In  the  present  case  of  a  planar  dipole,  the  length  of 
the  planar  dipole  is  selected  such  that  L/2w  =  33  [Ref.  9].  Figure  5.2  shows,  on  a 
semi— log  scale,  the  results  computed  by  the  developed  program  of  the  monostatic 
normalized  radar  cross  section  a/X  of  this  single  dipole  for  values  of  L/A  ranging 
from  0  to  1.4,  and  for  loads  Zl  =  0  and  Zl  =  oo.  The  case  of  Zl  =  oo  is  achieved  by 

setting  a  small  gap  (0.01  wavelength)  at  the  center  of  this  planar  dipole.  In  Figure 

2 

5.2,  the  corresponding  values  of  a/X    from  [Ref.  8,  pp.  115]  are  also  shown. 

This  comparison  leads  to  an  accurate  check  of  the  correct  calculation  of  the 
radar  cross  section  by  the  developed  RCS  program.  The  small  discrepanc;'  between 
the  computed  values  and  those  given  by  [Ref.  8]  is  due  to  the  error  incurred  in 
reading  values  off  the  curves  presented  in  [Ref.  8,  pp.  115]. 
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Figure  5.1   Vertical  Centered  Loaded  Planar  Dipole. 
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Figure  5.2  Normalized  RCS  a/ A*'  of  a  Center  Loaded  Planar  Dipole. 
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B.   SCATTERING  FROM  PLANAR  FRACTAL  TREES 

In  Chapter  4,  five  types  of  planar  fractal  trees  were  described.  Each  type  of 
these  trees  was  lying  in  the  y-z  plane.  The  developed  RCS  program  calculates  the 
broadside  radar  cross  section  a  of  each  tree  for  frequencies*  15  GHz-75  GHz,  for  the 
monostatic  case.  The  initial  dipole  that  is  being  reduced  by  the  reduction  factor  r 
has  physical  length  2  cm  and  only  the  branch  lengths  are  changed  for  each  tree 
depending  upon  the  reduction  factor  r.  As  the  frequency  of  the  plane  wave  changes, 
the  electrical  length  of  this  dipole  as  well  as  the  dipoles  composing  the  tree  will  also 
change.  For  the  frequency  range  of  15  GHz-75  GHz,  the  electrical  length  of  the 

initial  length  will  vary  from  one  to  four  wavelengths. 

2 
Figures  5.3—5.6    show  the  variations  of  a/X  ,  in  dB,  in  the  E— plane  and  the 

H— plane  for  the  case  of  a  fractal  tree  characterized  by  reduction  factor  r  =  0.53  and 

branch  angle  6  =  29.940.  The  frequency  of  the  incident  plane  wave  varied  from  15 

GHz  to  60  GHz  in  steps  of  15  GHz.  The  tree  is  composed  of  31  dipoles.  This  number 

is  small  as  the  reduction  factor  is  small  and  the  branch  lengths  are  reduced  to  very 

small  values  at  low  levels. 

9 

The  (t/A"  variations  in  the  E-Plane  are  in  a  range  of  10  dB  approximately.  In 

2 
the  H-Plane.  The  a/X    is  symmetric  about  the  90^  axis  and  is  smoother  at  lower 

2 
frequencies.  As  the  frequency  increases,  the  maximum  value  of  a/X    at  ^o  =  90^  and 

(^=00,  increases  from  2.73  dB  at  15  GHz  to  19.23  dB  at  60  GHz.  Figure  5.7  shows 

2 
the  variation  of  the  maximum  c/X    in  terms  of  frequency  increments  for  the  same 

structure.  This  variation  is  very  small  between  30  and  45  GHz. 
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Figure  5.3  Normalized  RCS  in  dB  at  F  =  15  GHz  of  a  Fractal  Tree 
with  r  =  0.53  and  Q  =  29.940. 
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Figure  5.4  Normalized  RCS  at  F  =  30  GHz  of  a  Fractal  Tree 
with  r  =  0.53  and  Q  =  29.940. 
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Figure  5.5  Normalized  RCS  at  F  =  45  GHz  of  a  Fractal  Tree 
with  r  =  0.53  and  Q  =  29.940. 
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Figure  5.6  Normalized  RCS  at  F  =  60  GHz  of  a  Planar  Fractal  Tree 
with  r  =  0.53  and  9  =  29.940. 
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Figure  5.7   Variation  of  maximum  cr/A    vs.  Frequency 

of  a  Planar  Fractal  Tree  with  r  =  0.53  and  6  =29.940. 
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Figures  5.8—5.12  show  the  variations  of  a/X"  ,  in  dB,  in  the  E  and  H  planes 
for  the  fractal  tree  characterized  by  a  reduction  factor  r  =  0.55  and  branch  angle  0 
=  46.430.  In  this  case,  the  fractal  tree  is  composed  of  63  dipoles  and  spreads  more 
than  the  previous  one,  and  the  lengths  of  the  branches  are  bigger  than  the  previous 
ones  (due  to  a  larger  reduction  factor  of  0.55  instead  of  0.53).  This  results  in  more 
variations  for  (j/A"  in  both  E  and  H  planes. 

This   model   is  investigated  at   a  frequency  range  of  15-75   GHz.    At   all 

2 

frequencies,  the  a/X    varies  in  a  range  of  10—12  dB  approximately.  In  the  H-Plane 

9 

the  variation  of  a/X"  is  symmetric  about  the  90^  axis.  The  maximum  value  of  radar 
cross  section  varies  from  2.78  db  at  15  GHz  to  22.63  dB  at  75  GHz.  Figure  5.13 

9 

shows  the  maximum  value  of  a/A",  corresponding  to  ^o  =  90^  and  v'o  =  0^  in  terms 
of  the  frequency  of  the  incident  plane  wave,  varying  from  15  GHz  to  75  GHz.  This 

9 

maximum  a/A"  increases  linearly  as  the  frequency  increases  except  the  range  of 
45-60  GHz  where  the  variation  is  verv  small. 
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Figure  5.8  Normalized  RCS  at  F  =  15  GHz  of  a  Planar  Fractal 
with  r  =  0.55  and  d  =  46.430. 
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Figure  5.9  Normalized  RCS  at  F  =  30  GHz  of  a  Planar  Fractal  Tree 
withr  =  0.55  and  0=  46.430. 
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Figure  5.10   Normalized  RCS  at  F  =  45  GHz  of  a  Planar  Fractal  Tree 
with  r  =  0.55  and  9  =  46.430. 
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Figure  5.1 1   Normalized  RCS  at  F  =  60  GHz  of  a  Planar  Fractal  Tree 
with  r  =  0.55  and  0  =^  46. 43^^. 
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Figure  5.12  Normalized  RCS  at  F  =  75  GHz  of  a  Planar  Fractal  Tree 
with  r  =  0.55  and  9  =  46.430. 
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Figure  5.13  Variation  of  Maximum  RCS  vs.  Frequency 
of  a  Planar  Fractal  Tree  with  r  =  0.55  and  Q  =  46.430. 
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Figures  5.14—5.18  show  the  variations  of  the  normalized  radar  cross  section 
cr/A"  of  a  planar  fractal  tree  characterized  by  reduction  factor  r  :=  0.60  and  branch 
angle  6  =  95.890.  xhis  model  is  composed  of  63  dipoles  and  has  more  spreading  than 
the  two  models  that  were  previously  investigated.  The  branches  have  larger  physical 

9 

lengths  than  the  ones  in  other  two  models.  At  15  GHz,  the  variation  of  a/X"  in  the 
E— Plane  is  very  low  and  similar  to  the  variation  in  the  other  two  models  at  the 
same  frequency.  In  the  H— Plane  this  variation  is  different  as  the  incident  plane 
wave  strikes  more  dipoles  from  —90^  to  90^. 

As  the  frequency  increases,  the  electrical  length  of  the  dipoles  is  also  increased, 
and  more  lobes  appear  at  30  GHz  and  higher  frequencies  in  both  the  E  and  the  H 
planes.  The  maximum  value  of  a/X  at  6o=  90^  and  tpo  =  0^  varies  from  -0.23  dB 
at  15  GHz  to  21.60  dB  at  75  GHz.  Figure  5.19  shows  the  variation  of  the  maximum 

radar  cross  section  in  terms  of  the  frequency  of  the  incident  plane  wave.  In  a  range 

2 
of  15—60  GHz  the  values  of  maximum  a/X    follow  a  straight  line  approximately.  In 

2 
the  range  of  60—75  GHz  the  variation  of  maximum  a/X    is  small. 

The  planar  fractal  trees  characterized  by  r  =  0.66,  $  =  145. 35^,  and  r  =  0.71,  0 

=  1800,  have  large  values  of  branching  angles.  The  large  values  of  reduction  factor  r 

generate  fractal  trees  whose  branches  have  large  physical  lengths  compared  with  the 

lengths  of  the  branches  of  the  previously  investigated  fractal  trees.  The  result  is  that 

the  electrical  lengths  of  these  branches  are  large  also  at  the  range  of  15—75  GHz,  and 

a  large  number  of  modes  is  required  to  investigate  the  last  two  types  of  fractal  trees. 

It  was  found  that  the  developed  RCS  program  is  not  able  to  calculate  the  radar 

cross  section  of  fractal  trees  that  are  characterized  by  large  values  of  reduction 

factor  r  and  branching  angle  $  due  to  memory  restrictions  and  other  numerical 

problems. 
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For  this  reason  the  numerical  results  for  ajX"  are  not  presented  in  this  thesis  for 
these  two  cases. 

It  is  seen  by  comparing  Figures  5.7,  5.13,  and  5.19  that  there  is  a  frequency 
range  over  which  the  variation  of  maximum  RCS  is  rather  small.  Furthermore,  this 
range  of  frequencies  shifts  to  higher  frequencies  as  the  tree  structure  is  spread  from  r 
=  0.53  to  r  =  0.60.  The  maximum  RCS  of  each  of  these  trees  varies  in  a  range  of  3 
dB  approximately  at  the  same  frequency. 

Scattering  from  an  actual  tree  will  not  exactly  follow  all  these  patterns,  but  it 
is  felt  that  the  trends  would  generally  remain  the  same.  This  is  specially  true  for  the 
frequency  behavior. 
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Figure  5.14  Normalized  RCS  at  F  =  15  GHz  of  a  Planar  Fractal  Tree 
with  r  =  0.60  and  0  =  95.890. 
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Figure  5.15  Normalized  RCS  at  F  =  30  GHz  of  a  Planar  Fractal  Tree 
with  r  =  0.60  and  $  =  95.890. 
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Figure  5.16   Normalized  RCS  at  F  =  45  GHz  of  a  Planar  Fractal  Tree 
with  r  =  0.60  and  0  =  95. 89^. 
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Figure  5.17   Normalized  RCS  at  F  =  60  GIIz  of  a  Planar  Fractal  Tree 
with  r  =  0.60  and  0  =  95.890. 


61 


20 

h 

2 

10 
0 

T\  ^m 

A 

^ 

w        f^ 

b 

-10 

M 

-20 

45 


90 

Go 


F    =    75    GHz 
E-Plane 


135 


180 


i^; 


J5U 

25 

20 
15 

A 

\ 

A 

^ 

10 

^\ 

\\r 

b 

5 
0 

V 

11 

F    =    75    GHz 
H-Plane 


-90 


-45 


45 


90 


(pO 
Figure  5.18  Normalized  RCS  at  F  =  75  GHZ  of  a  Planar  Fractal  Tree 
with  r  =  0.60  and  $  =  95.890. 
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Figure  5.19  Variation  of  Maximum  RCS  vs.  Frequency 
of  a  Planar  Fractal  Tree  with  r  =  0.60  and  ^  =  95.890. 
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VI.    CONCLUSIONS  AND  RECOMMENDATIONS 

The  radar  cross  section  of  natural  trees  is  a  complicated  problem  of 
electromagnetic  scattering.  The  real  trees  are  made  up  of  long,  intersecting,  and 
lossy  objects.  The  geometry  of  these  objects  is  not  easy  to  set  up  in  a  three 
dimensional  space. 

In  this  thesis,  the  real  fractal  trees  were  approximated  by  planar  fractal  trees. 
The  radar  cross  section  of  these  trees  was  calculated  by  developing  a  Fortran 
program.  The  planar  fractal  trees  that  were  used  for  this  investigation  were 
symmetric  planar  structures  composed  of  perfectly  conducting  and  non— intersecting 
planar  dipoles.  The  geometry  of  these  structures  was  generated  using  fractal  theory. 

A.      CONCLUSIONS 

The  radar  cross  section  of  the  planar  fractal  trees  was  calculated  using  the 
moment  method  theory.  For  a  given  planar  structure  composed  of  planar  strips  and 
illuminated  by  a  plane  wave,  the  program  generates  the  geometry  of  the  PWS 
modes,  calculates  the  impedance  between  any  two  of  them,  and  fills  the  impedance 
matrix  [Z].  The  voltage  on  each  PWS  mode,  due  to  the  incident  electric  field,  is  also 
calculated,  and  the  voltage  column  vector  [V]  is  generated. 

The  RCS  program  uses  the  moment  method  theory  to  calculate  the  current 
distribution  on  each  PWS  mode.  The  radiation  equations  of  electromagnetic  theory 
are  used  to  determine  the  scattered  electric  field  due  to  the  current  induced  on  each 
PWS  mode  of  the  structure.  The  knowledge  of  both  the  given  incident 
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electric  and  calculated  scattered  electric  fields  leads  to  the  calculation  of  the  radar 
cross  section   [ajX  ),  for  the  given  angles   ^o,V^o  of  the  incident   electric  field. 

The  calculated  radar  cross  section  of  a  single  centered  loaded  vertical  dipole 
was  compared  with  the  values  given  in  [Ref.  8,  pp.  115],  for  a  similar  dipole. 
Although  the  investigation  was  for  a  planar  dipole  and  in  [Ref.  8]  the  dipole  was 
cylindrical,  the  discrepancy  of  the  results  was  very  small. 

In  this  thesis,  five  models  of  planar  fractal  trees  were  investigated.  The 
geometry  of  these  models  was  generated  by  a  given  program  which  was  modified  to 
generate  the  input  data  for  the  developed  RCS  program. 

The  broadside  radar  cross  section,  for  the  monostatic  case,  was  calculated  for 
three  of  these  planar  fractal  trees,  for  a  frequency  range  of  15  GHz  to  75  GHz  in 
steps  of  15  GHz.  This  investigation  showed  that  the  radar  cross  section  varied  in  a 

9 

range  of  10  dB  in  the  E— Plane.  In  the  H— Plane,  the  variation  of  ajX"  was 
symmetric  about  the  90^  axis  and  smooth  at  low  frequencies  and  small  branching 
angles.  For  higher  branching  angles  and  reduction  factors,  more  variations  of  cr/A" 
with  the  frequency  were  seen  in  both  the  E  and  H  planes. 

The  maximum  RCS  at  ^o  =  90^,  ip^  =  QO  showed  a  small  variation  over  a 
frequency  range.  This  range  was  different  on  each  investigated  tree.  The  variation  of 
the  maximum  RCS  followed  a  straight  line  at  the  other  frequencies.  In  each  tree  and 
at  the  same  frequency,  the  maximum  RCS  varied  in  the  range  of  3  dB 
approximately. 

It  was  found  that  the  developed  RCS  program  was  not  able  to  calculate  the 
radar  cross  section  of  fractal  trees  that  were  characterized  by  large  values  of 
reduction  factor  and  branching  angle  due  to  memory  restrictions  and  other 
numerical  problems. 
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The  scattering  from  an  actual  tree  will  not  exactly  follow  the  patterns  that 
were  described  in  Chapter  5,  but  it  is  felt  that  the  trends  would  generally  remain 
the  same.  This  is  specially  true  for  the  frequency  behavior. 

B.      RECOMMENDATIONS 

One  recommendation  is  to  investigate  the  radar  cross  section  of  planar  fractal 
trees  characterized  by  values  other  than  those  used  in  this  thesis.  Especially,  a 
limited  set  of  reduction  factor  and  branching  angle  has  to  be  established  for  the 
same  physical  length  of  the  initial  dipole. 

Another  recommendation  is  to  investigate  the  radar  cross  section  of  planar 
fractal  trees  characterized  by  the  same  values  of  reduction  factor  and  branching 
angles  as  those  that  were  used  in  this  thesis  but  with  different  physical  and 
electrical  dimensions  of  the  fractal  trees. 

In  this  thesis,  the  branches  of  the  fractal  trees  were  assumed  to  be  perfectly 
conducting  planar  strips.  The  trees  were  considered  without  leaves.  The  radar  cross 
section  of  fractal  trees  composed  of  lossy  planar  strips  and  with  leaves  should  be 
investigated. 

Finally,  the  generation  of  a  three  dimensional  fractal  tree  and  its  radar  cross 
section  may  be  investigated. 
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APPENDLX  A 
PROGR„\M  RCS 


C 

C  THIS  PROGRAM  CALCULATES  THE  RADAR  CROSS  SECTION  OF  A 

C  PLANAR  STRUCTURE  COMPOSED  OF  ARBITRARILY  ORIENTED 

C  PLANAR  THIN  DIPOLES.  THE  DIPOLES  ARE  NOT  INTERSECTING. 

C  THE  DIPOLES  ARE  PERFECT  CONDUCTORS. 

C 

c 

C  GEOMETRY  IN  Y-Z  PLANE 

C 

C  INPUT  DATA: 

C 

C  F  =  FREQUENCY  IN  GHZ 

C  N\V  -  NUMBER  OF  DIPOLES 

C  A.B  =  COORDINATES  OF  INCIDENT  ELECTRIC  FIELD  ON  Y  AND  Z 

C  AXIS  RESPECTIVELY. 

C  L\V   =  HALF  LENGTH  OF  EACH  DIPOLE  IN  CM 

C  WW  =   HALF  WIDTH  OF  EACH  DIPOLE  IN  CM 

C  SW.  TW  =  COORDINATES  OF  THE  CENTER  OF  EACH  DIPOLE 

C  ALONG  Y  AND  Z  AXIS  RESPECTRTLY 

C  PSIW  =  ANGLE  IN  DEGREES  BETWEEN  THE  Z  AXIS  (REFERENCE) 

C  AND  DIRECTION  OF  CURRENT  FLOW  ON  EACH   DIPOLE 

C  ( POSITIVE  ANGLES  ARE  MEASURED  CLOCKWISE ) 

C  NW  =  NUMBER  OF  SEGMENTS  THAT  EACH  DIPOLE  IS 

C  SUBDRIDED 

C  THITAO.  PHIO  =  ANGLES  IN  DEGREES  OF  INCIDENT  ELECTRIC 

C  FIELD. 

C 

C  OUTPUT  DATA: 

C 

C  NORMALIZED  RCS  ia/X')  IN  dB 

C 

REAL   LS(1:230).S(1:230).T(1:230).PSI(1:230) 

REAL   LEN.PG.SU.TU.WID.A.B.MAG(1:230) 

INTEGER  I.NMAXTNDXf230).NT.NP.L.M.N.K.NW 

COMPLEX  Z(l:230.1:230).V(l:230j.CUR(T:230j 

COMMON/  RCl/  NMAX 

COMMON/  RC2/  LSAVLS.T.PSI 

COMMON/  RC4/  CUR 

OPEN  (UNIT=1.FILE=TNPUT1'.F0RM='F0RMATTED') 

0PEN(UNIT=2.FILE='0UTPUTr.F0RM='UNF0RMATTED'j 

OPEN  (UNIT=3.FILE=TNPUT'.F0RM='F0RMATTED' ) 

0PEN(UNIT=4.FILE='0UTPUT'.F0RM=T0RMATTED') 
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READ(1,*)  F,NW,A,B 

PI  =  4.*  ATAN(1.) 

KO  =  PI*F/15. 

RAD  =  PI/180. 

CALL  ZMATR  (F,NW,Z) 

NP-230 

NT=NMAX 

CALL  CLUDCP  (Z,NT,NP,INDX) 
22       READ  (3,*,END=11)  THITAO,PHIO 

CALL  VOLT(F,THITA0,PHI0,A,B,V) 

CALL  CLUBSB  (Z,NT,NP,INDX,V) 

DO  54  K=1,NT 

CUR(K)=V(K) 
54        CONTINUE 

CALL  RADIAT  (F,THITAO,PHIO,NMAX,A,B,RCS) 

WRITE(4,*)  THITAO,PHIO,RCS 

PRINT*,RCS 

GOTO  22 
11        CL0SE(3) 

STOP 

END 
C 

C  SUBROUTINE  TO  COMPUTE  THE  MUTUAL  IMPEDANCE  BETWEEN 

C  THE  PWS  MODES  OF  N  ARBITRARILY  ORIENTED  DIPOLES. 

C 

C  GEOMETRY  IN  THE  Y-Z  PLANE 

C 

SUBROUTINE  ZMATR  (F,NW,Z) 

REALF,S(1:230),T(1:230),LG,SG,TG,WI(1:230),PSI(1:230) 

REALPSIG,PSIW(1:70),LW(1:70),H,W,DH,DW,WIDTH,LENG 

REALH1,H2,W1,W2,S1,S2,T1,T2,PSI1,PSI2,LS(1:230) 

REALWW(1:70),SW(1:70),TW(1:70),SS(1:20),TS(1:20) 

REAL  PG,WID,LEN,SU,TU,A,B 

INTEGER  TEMP,P,L,NS(1:70),NMAX,I,N,NW 

INTEGER  FUN,G,J,GB,GR,K 

COMPLEX  Z12,Z(1:230,1:230) 

COMMON/  RC2/  LS,WI,S,T,PSI 

COMMON/  RCl/  NMAX 

COMMON/  JOHN/  LG,NG,SG,TG,PSIG,LENG 

COMMON/  ZPAR/  H,W,DW,DH 

COMMON/  ZINC/  H1,H2,W1,W2,S1,S2,T1,T2,PSI1,PSI2 

DO  101  I  =  1,  NW 
READ(1,*)LW(I),WW(I),SW(I),TW(I),PSIW(I),NS(I) 
101       CONTINUE 

CLOSE(l) 

NMAX  =  0 

DO  75  I=1,NW 
NMAX  =  NMAX  +  (NS  (I)  -  1) 
75       CONTINUE 

P=l 

GB=0 
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DO  50  1=1, NW 

LENG  -  LW(I) 

LG=2='LW(I)/FL0AT(NS(I)) 

WIDTH=W\V(I) 

NG=NS(I) 

SG=SW(I) 

TG=TW(I) 

PSIG=PSIW(I) 

CALL  GEOM(SS,TS) 

GR=NS(I)+GB-1 

M  =  P 

DO60  J=1,NS(I)-1 
PSI(M)=PSIG 
LS(M)=LG 
WI(M)=\VIDTH 
S(M)=SS(J) 
T(M)=TS(J) 
PRINT^S(M),T(M) 
M  =  M  +  1 
60  CONTINUE 

GB=GR 

P=GR+1 
50        CONTINUE 

TEMP=NS(1)-1 
L=l 

G=NMAX 
DO80M  =  l.G 

IF(M.LE.TEMP)GOTO  85 

L=L+1 

TEMP=TEMP+NS(L)-1 
85  CONTINUE 

K=TEMP 

DO  90  N=M.G 

IF(N.LE.K)THEN 

H=LS(M) 

W=WI(M) 

DW=0 

DH=ABS(M-N)*H 

CALL  ZSDIP(F,Z12) 

Z(M.N)=Z12 

PRINT^Z(M,N) 

WRITE(2)  Z(M,N) 

ELSE 

H1=LS(M) 

H2=LS(N) 

W1=WI(M) 

W2=WI(N) 

S1=S(M) 

S2=S(N) 

T1=T(M) 
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T2=T(N) 

PSI1=PSI(M) 

PSI2-PSI(N) 

CALL  ZPSUR  (F,Z12) 

Z(M,N)=Z12 

PRINT*,Z(M.N) 

WRITE(2)  Z(M,N) 

ENDIF 
90  CONTINUE 

80       CONTINUE 
DO  24  M=2,G 

DO  26  N=1,M-1 

Z(M,N)=Z(N,M) 

WRITE(2)  Z(M,N) 
26  CONTINUE 

24        CONTINUE 
RETURN 
END 
C 

C  SUBROUTINE  TO  COMPUTE  THE  MUTUAL/SELF  IMPEDANCE 

C  BETWEEN  TWO  DIPOLES.  THE  DIPOLES  ARE  ASSUMED  TO  BE 

C  COPLANAR,  IDENTICAL  AND  PARALLEL.  THE  DIPOLE  TO 

C  DIPOLE  IMPEDANCE  IS  COMPUTED  AS  THE  SUM  OF  FOUR 

C  MONOPOLE  TO  MONOPOLE  IMPEDANCES. 

C 

C  INPUT  PARAMETERS: 

C 

C  F  =  Frequency  of  operation  (GHz) 

C  H  =  Half  height  of  the  dipole  (cm) 

C  W  =  Half  width  of  the  dipole  (cm) 

C         DH  =  Longitudinal  distance  between  the  two  dipoles  (cm) 
C         DW  =  Transverse  distance  between  the  two  dipoles  (cm) 


C 


SUBROUTINE  ZSDIP  (F,Z12) 

REAL  H,  W,  DW,  DH,  F 

COMPLEX  Z12,  ZT 

COMMON/  ZPAR/  H,W,DW,DH 

ZT  =  (0.,0.) 

Z12  =  (0.,0.) 

CALL  ZSMONP  (F,  H,  W,  DW,  DH,  0,  1,  0,  1,  ZT) 

Z12  =  Z12  +  ZT 

CALL  ZSMONP  (F,  H,  W,  DW,  DH  +  H,  0,  1,  1,  0,  ZT) 

Z12  =  Z12  +  ZT 

CALL  ZSMONP  (F,  H,  W,  DW,  DH  -  H,  1,  0,  0,  1,  ZT) 

Z12  =  Z12  +  ZT 

CALL  ZSMONP  (F,  H,  W,  DW,  DH,  1,  0,  1,  0,  ZT) 

Z12  =  Z12  +  ZT 

RETURN 

END 
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C  SUBROUTINE  TO  COMPUTE  THE  SELF/MUTUAL  IMPEDANCE 

C  BETWEEN  TWO  IDENTICAL,  COPLANAR  MONOPOLES.  THE 

C  CURRENT  IS  ASSUMED  TO  BE  CONSTANT  IN  THE 

C  TRANSVERSE  DIRECTION. 

C 

C  REE:    R.  JANASWAMY,  A  SIMPLIFIED  EXPRESSION  FOR  THE 

C  SELF/MUTUAL  IMPEDANCE  BETWEEN  TWO  COPLANAR 

C  AND  PARALLEL  MONOPOLES,  IEEE  T-AP,  AP-35, 

C  No.  10,  pp.  1174-1176,  October  1987. 

C 

C  INPUT  PARAMETERS: 

C 

C    F  =  FREQUENCY  IN  GHz 

C    H  =  LENGTH  OF  EACH  MONOPOLE  (cm) 

C    W  =  WIDTH  OF  EACH  MONOPOLE  (cm) 

C    D  =  CENTER  TO  CENTER  SPACING  BETWEEN  THE  TWO 

C    MONOPOLES  IN  THE  DIRECTION  TRANSVERSE  TO  THE 

C    CURRENT  FLOW  (cm) 

C    HH  =  CENTER  TO  CENTER  SPACING  BETWEEN  THE  TWO 

C    MONOPOLES  IN  THE  DIRECTION  OF  CURRENT  FLOW  (cm) 

C 

C    111  =  TERMINAL  CURRENT  OF  END  1  OF  MONOPOLE  1. 

C    121  =  TERMINAL  CURRENT  OF  END  2  OF  MONOPOLE  1. 

C    112  =  TERMINAL  CURRENT  OF  END  1  OF  MONOPOLE  2. 

C    122  =  TERMINAL  CURRENT  OF  END  2  OF  MONOPOLE  2. 

C 

C         NOTE:  111,  121, 112,  122  can  assume  values  onlvO  or  1. 

C  ICODE  =  0.  IF   D   .LE.  4W 

C  ICODE  =  1,  OTHERWISE 

C  With  ICODE  =  0,  the  expression  provided  in  the  above  paper  is  used. 

C  With  ICODE  =  1,  a  modified  form  of  the  expression  provided  in  the 

C         above  paper  is  used.  (cf.  notes) 

C 

C  OUTPUT  PARAMETERS: 

C 

C  Z12  =  COMPLEX  IMPEDANCE  BETWEEN  THE  TWO  SURFACE 

C  MONOPOLES. 

C 

SUBROUTINE  ZSMONP  (F,  H,  W,  D,  HH,  111,  121,  112,  122,  Z12) 

REAL  F,  H,  W,  D,  HH,  A,  B,  PI,  KO,  V,  UB,  UBP,  UA,  UAP 

REAL  KW,  KH.  KD,  RC,  FR,  EI,  II,  12,  13,  14,  SI,  CI 

REAL  UABP.  ZR,  ZI,  X,  AA  (1),  BB  (1),  SQXV,  UAB,  TINY 

INTEGER  111,  121,  112,  122,  M,  N,  NX,  KI,  ICODE 

COMPLEX  Z12,  AC  (-1:1,  -1:1),  El,  E2,  Zl,  J,  CMN,  E3,  EI,  FAC 

EXTERNAL  FR.  FI 

COMMON  /PARAM/  N,  V,  KD,  A,  B,  ICODE 

EI(X)  =  CI(ABS(X))-J*SI(X) 

SQXV  (X)  =  SQRT  (X  =*  X  +  V  *  V) 

TINY  =  l.E-6 

PI  =  4.  *  ATAN  (1.) 

NX  =  1 


KO  =  PI  *  F  /  15. 

KW  =  KO  *  W 

KH  =  KO  *  H 

KD  =  KO  *  D 

ICODE  =  0 

IF  (KD  .GT.  4.  *  KW)  ICODE  =  1 

KI  =  10 
C 

C         Note  that  with  this  choice  of  KI  accurate  results  are  found  in  the 
C  region  0  <  D  <  4  W,  and  for  D  >  20  W.  In  between  these  two  regions, 

C  accurate  results  are  found  with  KI  =  20.  Hence  the  above  choice  is 

C         good  only  if  this  code  is  used  for  values  of  D  satisfying  the  above 
C  inequalities. 


C 


FAC  =  CMPLX(1.,0.) 

A  =  KD-2.  *KW 

B  =  KD  +  2.  *  KW 

IF  (ICODE  .EQ.  1)  GO  TO  2 

AA  (1)  =  A 

BB  (1)  =  B 

GO  TO  3 

AA  (1)  =  -2.  *  KW 

BB  (1)  =  2.  *  KW 

II  =  (121  *C0S(KH)-I11) 

12=    111-121  *  COS  (KH) 

13=    121  -111  *  COS   KH 

14  =  (121  -111  *  COS  (KH) 

J  =  CMPLX  (0.,1.) 

Zl  =  CEXP  (J  *  KH) 

AC  (-1,-1)  =  12 +  11  /Zl 

AC  (-1,  1)  =  12  +  II  *  Zl 

AC    1,-1)  =  13-14*  Zl 

AC    1,  1)  =  13-14/  Zl 

AC   0,  -1)  =  -  (AC  (-1,  -1)  *  Zl  +  AC  (1,  -1)  /  Zl) 

AC  (0,  1)  =  -  (AC  (1,  1)  *  Zl  +  AC  (-1,  1)  /  Zl) 

RC  =  15.  /  (2.  *  SIN  (KH)  *  KW)  **  2 

Z12  =  (0.,0.) 

DO  1  M  =  -1,  1 

DO  1  N  =  -l,  1,2 

V  =  KO  *  (HH  +  M  *  H) 

IF  (ICODE  .EQ.  0)  GO  TO  4 

FAC  =  CEXP  (J  *  N  *  V) 

CMN  =  CMPLX  (0.,0.) 

GO  TO  5 

UA  =  SQXV  (A)  4-  N  *  V 

UAP  =  UA-2.  *N*  V 

UB  =  SQXV  (B)  +  N  *  V 

UBP  =  UB-2.  *N*  V 

UAB  =  SQXV  (KD)  +  N  *  V 

UABP  =  UAB-2.  *N*  V 

El  =  EI  (UB) 
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E2  =  (0.,0.) 

IF  (ABS  (A)  .GT.  TINY)  E2  =  EI  (UA) 

E3-  (0..0.) 

IF  (ABS  (KD)  .GT.  TINY)  E3  =  EI  (UAB) 

CMN  =  0.5  *  (B  *  B  *  El  +  A  *  A  *  E2  -  2.  *  KD  *  KD  *  E3  + 
k  CEXP  (-J  *  UB)  *  (1.  +  J  *  UBP)  +  CEXP  (-J  *  UA)  * 

k  (1.  +  J  *  UAP)  -  2.  *  CEXP  (-  J  *  UAB)  *  (1.  +  J  *  UABP)) 

5        CALL  HABER  (NX.  AA,  BB,  FR,  KI,  ZR) 

CALL  HABER  (NX.  AA.  BB,  FI,  KI,  ZI) 

Z12  =  Z12  +  AC  (M,  N)  *  (CEXP  (J  *  N  *  V)  *  CMN  +  (ZR  +  J  *  ZI) 
k  *  FAC) 

1        CONTINUE 

Z12  =  Z12  =^  RC 

RETURN 

END 
C 

REAL  FUNCTION  FR  (XI,  NX) 

INTEGER  N.  NX,  CODE 

REAL  X.  V,  Tl.  KD.  A,  B.  XI  (NX),  TKW,  T2,  CI.  SI,  TINY 

COMPLEX  EI.  J 

COMMON  /PARAM/  N,  V,  KD.  A,  B,  CODE 

EI(X)  =  CI(ABS(X))-J*SI  (X) 

TINY  =  LE-6 

X  =  XI  (1) 

J  =  CMPLX(O..L) 

IF  (CODE.EQ.  1)  GOTO  1 

Tl  =  SQRT  (X  *  X  +  V  *  V) 

FR  =  C0S(T1) 

IF  (ABS  (V)  .GT.  TINY)  FR  =  FR  *  (1.  -  N  *  V  /  Tl) 

IF  (X  .LE.  KD)  FR  =  FR*  A 

IF  (X  .GT.  KD)  FR  =  -FR  *  B 

GO  TO  2 

1  TKW=   0.5*(B-A) 

T2  =  SQRT  ((KD  +  X)  *  (KD  +  X)  +  V  *  V) 
FR  =  REAL  (EI  (T2  +  N  ^  V)) 
FR  =  FR*  (TKW -ABS  (X)) 

2  END 


REAL  FUNCTION  FI  (XI,  NX) 

INTEGER  N,  NX,  CODE 

REAL  X,  V,  Tl,  KD,  A,  B,  XI  (NX),  TKW,  T2,  CI,  SI,  TINY 

COMPLEX  EI,  J 

COMMON  /PARAM/  N,  V,  KD,  A,  B,  CODE 

EI(X)  =  CI  (ABS  (X)) -J*  SI  (X) 

X  =  XI(1) 

TINY=  l.E-6 

J  =  CMPLX(0.,  1.) 

Tl  =  SQRT  (X  *  X  +  V  *  V) 

IF  (CODE.EQ.  1)  GOTO  1 
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FI  =  -SIN  (Tl) 

IF  (ABS  (V)  .GT.  TINY)  FI  =  FI  *  (1.  -  N  *  V  /  Tl) 

IF  X  .LE.  KD)  FI  =  FI  *  A 

IF  (X  .GT.  KD)  FI  =  -FI  *  B 

GO  TO  2 

1  TKW=  0.5*  (B-A) 

T2  =  SQRT  ((KD  +  X)  *  (KD  +  X)  +  V  *  V) 
FI  =  AIMAG  (EI  (T2  +  N  *  V)) 
FI  =  FI*  (TKW-ABS(X)) 

2  END 
C 

C  SUBROUTINE  TO  COMPUTE  THE  COORDINATES  OF  THE  PWS 

C  MODES  OF  EACH  DIPOLE  IN  THE  Y-Z  PLANE 

C 

SUBROUTINE  GEOM  (S,T) 

REAL  THI,H1,H2,Y,Z,SG,TG,LG,PSIG 

REAL  P,Q,S(1:20),T(1:20),LENG 

INTEGER  M,NG,K 

COMMON/  JOHN/  LG,NG,SG,TG,PSIG,LENG 

COSD(X)  =  COS(X*RAD) 

SIND(X)  =  SIN(X*RAD) 

PI=4.*ATAN(1.) 

RAD=PI/180. 

THI=90.-PSIG 

H1=C0SD(THI) 

H2=SIND(THI) 

Y=SG-(LENG*H1) 

Z=TG-(LENG'H2) 

P=LG*H1 

Q=LG*H2 

S(1)=Y+P 

T(1)=Z+Q 

DO  225  K=2,NG-1 

S(K)=S(1)+(K-1)*P 
T(K)-T(1)+(K-1)*Q 
225       CONTINUE 

RETURN 

END 


SUBROUTINE  ZPSUR  (F,Z12) 

REAL  HI,  Wl,  H2,  W2,  PSI,  F,  YSTAR,  ZSTAR,  KO,  AI  (2),  BI  (2) 

REAL  C  (3),  D  (3),  RT,  ZT,  COT,  CSEC,  KOS,  X,  PI,  PSIl,  PSI2 

REAL  SI,  Tl,  S2,  T2,  RINTG,  IMINTG,  RESULT,  SIND,  COSD,SI,CI 

INTEGER  NX,  KI 

COMPLEX  Z12,  J 

EXTERNAL  RINTG,  IMINTG 

COMMON/  ZINC/  H1,H2,W1,W2,S1,S2,T1,T2,PSI1,PSI2 

COMMON  /PARAM2/  C,  D,  RT,  ZT,  CSEC,  COT,  KOS,  J,  KO 

SIND  (X)  =  SIN  (X  *  PI  /  180.) 


74 


COSD  (X)  =  COS  (X  ^  PI  /  180.) 

DATA  AL  BI/2  =--!..  2^  1./ 

PI  =  4.  =- ATAN  (1.) 

PSI  =  PSI2-PSI1 

IF  (ABS  (PSI).LE.0.4)  PSI=  SIGN(1.,PSI)*0.4 

J  =  CMPLX  (O.J.) 

NX  =  2 

KI  =  3 

KO  =  PI  ^  F  /  15. 

C(l)  =  1.  /SIN(K0*H1) 

C  (3)  =  C  (1) 

C(2)  =  -2.  =^COS(KO*H1)  *C(1) 

D  (1)  =  1.  /  SIN  (KO  "-  H2) 

D(3)  =  D(1) 

D  (2)  =  -2.  *  COS  (KO  "^  H2)  *  D  (1) 

CSEC  =  1.  /  SIND  (PSI) 

KOS  =  COSD  (PSI) 

COT  =  KOS  =-  CSEC 

YSTAR  =  (S2-S1)  '  COSD  (PSIl)  -  (T2-T1)  *  SIND  (PSIl) 

ZSTAR  =  (S2-S1)  '^  SIND  (PSIl)  +  (T2-T1)  '  COSD  (PSIl) 

RT  -  YSTAR  -  CSEC  -  H2 

ZT  =  YSTAR  *  COT  -  HI  -  ZSTAR 

CALL  HABER  (NX.  AL  BI,  RINTG.  KI,  RESULT) 

Z12  =  RESULT 

CALL  HABER  (NX,  AL  BI.  IMINTG.  KI,  RESULT) 

Zr2  =  Z12  +  J  ="  RESULT 

Z12  =  -3.75  -^  Z12 

RETURN 

END 

REAL  FUNCTION  SI  (XI) 

REAL  XI 

REAL  AF  (4).  BE  (4).  AG  (4).  BG  (4).  X.  X2,  X4.  X6.  X8.  EX.  GX, 

PL  SGN 
DATA   AF  /  .38.027264.  265.187033.  335.67732.  38.102495  / 
DATA   BE  /  40.0214.33.  322.624911.  570.23628,  157.105423  / 
DATA   AG'/  42.242855.  .302.757865.  352.018498.  21.821899  / 
DATA   BG  /  48.196927,  482.485984.  1114.978885,  449.690326  / 
SI  =  0. 

IF  (XI  .EQ.  0.)  RETURN 
SGN  =  +1. 

IF  (XI  .LT.  0.)  SGN  =  -1. 
X  =  ABS(XI) 
X2  =  X  -  X 
IF  (X  .GE.  20.)  THEN 
FX  =  \.  /X'  (1.-2.  /X2) 
GX  =  1.  /X2''  (1.-6.  /X2) 
GOTO  1 
END  IF 
X4  =  X2  *  X2 
X6  =  X4  *  X2 
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X8  =  X4  *  X4 

IF  (X  .LT.  1.)  THEN 

SI  =  X*(l.  -  X2  /  18.  +  X4  /  600.  -  X6  /  35280.  +  X8  /  3265920.) 

SI  =  SI  ^  SGN 

ELSE 

FX  =  (1.  +  AF  (1)  /  X2  +  AF  (2)  /  X4  +  AF  (3)  /  X6  +  AF  (4)  / 
k  X8) 

FX  =  FX  /  (X  *  (1.  +  BF  (1)  /  X2  +  BF  (2)  /  X4  +  BF  (3)  /  X6  + 
k  BF  (4)  /  X8)) 

GX  =  (1.  +  AG  (1)  /  X2  +  AG  (2)  /  X4  +  AG  (3)  /  X6  +  AG  (4)  / 
k  X8) 

GX  =  GX  /  (X2  *  (1.  +  BG  (1)  /  X2  +  BG  (2)  /  X4  +  BG  (3)  /  X6  + 
k  BG(4)  /  X8)) 

1        PI  -  4.  *  ATAN  (1.) 

X  =  X-AINT(X/(2.  *PI))  *2.  *PI 

SI  =  SGN  *  (PI  /  2.  -  FX  *  COS  (X)  -  GX  *  SIN  (X)) 

END  IF 

END 


REAL  FUNCTION  CI  (X) 

REAL  AF  (4),  BF  (4),  AG  (4),  BG  (4),  X,  X2,  X4,  X6,  X8,  FX,  GX 

REAL  PI 

DATA   AF  /  38.027264,  265.187033,  335.67732,  38.102495  / 

DATA   BF  /  40.021433,  322.624911,  570.23628,  157.105423  / 

DATA   AG  /  42.242855,  302.757865,  352.018498,  21.821899  / 

DATA   BG  /  48.196927,  482.485984,  1114.978885,  449.690326  / 

IF  (X  .LE.  0.)THEN 

PRINT  *,  'Invalid  argument  for  CI  (x)',  '  x  =  ',  x 

RETURN 

ELSE 

X2  =  X  *  X 

X4  =  X2  *  X2 

IF  (X  .GE.  20.)  THEN 

FX  =  1.  /X*  (1.-2.  /X2) 

GX  =  1.  /X2*(l--6.  /X2) 

GOTO  1 

END  IF 

X6  =  X4  *  X2 

X8  =  X4  *  X4 

IF  (X  .LT.  1.)  THEN 

CI  =  0.57721566  +  ALOG  (X)  -  X2  *  (0.25  -  X2  /  96.  +  X4  /  4320. 

-X6/ 322560.) 
ELSE 
FX  =  (1.  +  AF  (1)  /  X2  +  AF  (2)  /  X4  +  AF  (3)  /  X6  +  AF  (4)  / 

X8) 
FX  =  FX  /  (X  *  (1.  +  BF  (1)  /  X2  +  BF  (2)  /  X4  +  BF  (3)  /  X6  + 

BF  (4)  /  X8)) 
GX  =  (1.  +  AG  (1)  /  X2  +  AG  (2)  /  X4  +  AG  (3)  /  X6  +  AG  (4)  / 

X8) 
GX  =  GX  /  (X2  *  (1.  +  BG  (1)  /  X2  +  BG  (2)  /  X4  +  BG  (3)  /  X6  + 
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k  BG(4)  /  X8)) 

1        PI  =  4.  *  ATAN  (1.) 

X  =  X  -  AINT  (X  /  (2.  *  PI))  *  2.  *  PI 

CI  =  FX  *  SIN  (X)  -  GX  *  COS  (X) 

END  IF 

END  IF 

END 
C 

REAL  FUNCTION  RINTG  (X,  NX) 

COMPLEX  J,  EI,  TERMl,  TERM2,  INTGR 

INTEGER  NX.  M.  N.  P.  Q 

REAL   '  4  PSIl,  PSr2,  Si.  Tl,  S2,  T2 

REAL  =^  4  C  (3),  D  (3),  R,  Z,  RT,  ZT,  U,  V,  CSEC,  COT,  KOS.  Y 

REAL  *  4  X  (NX),  HI.  H2,  Wl,  W2,  KO,  PZQR,  SI,  CI 
k  .TT,  ZM.  RN.  RMN,  TESC,  TESD,  PI,  PZQRP 

COMMON/  ZINC/  H1,H2,W1,W2,S1,S2,T1,T2,PSI1,PSI2 

COMMON  /PARAM2/  C,  D.  RT,  ZT,  CSEC,  COT,  KOS,  J,  KO 

EI  (Y)  =  CI  (ABS(Y))-J=^SI(Y) 

U  =  X(1) 

V  =  X  (2) 

INTGR  =  (0..0.) 

TESC  =  ABS  (C  (2)  /  C  (1)) 

TESD  =  ABS(D  (2)  /  D  (1)) 

TT  =  ALOG  (ABS  ((l.+KOS)/(l.-KOS))) 

PI  =  4.  ="  ATAN  (1.) 

DO  3  M  =  1,3 

IF  (M  .EQ.  2  .AND.  TESC  .LE.  l.E-6)  GO  TO  3 

Z  =  (M-1)  *  HI 

ZM  =  -U  *  Wl  *  COT  +  V  *  W2  *  CSEC  +  ZT  +  Z 

TERM2  =  (0.,0.) 

D0  2N=  1,3 

IF  (N  .EQ.  2  .AND.  TESD  .LE.  l.E-6)  GO  TO  2 

R  =  (N-1)  "^  H2 

RN  =  -U  *  Wl  *  CSEC  +  V  *  W2  *  COT  +  RT  +  R 

IF  (ABS  (KO  *  RN)  .LE.  0.5E-1)  THEN 

PZQRP  =  KO  *  ABS  (ZM)  -  AINT  (KO  *  ABS  (ZM)  /  (2.  *  PI))  * 
k  2.  *  PI 

TERMl  =  CEXP  (-J  *  PZQRP)  *  TT 

GO  TO  4 

END  IF 

IF  (ABS  (KO  *  ZM)  .LE.  0.5E-1)  THEN 

PZQRP  =  KO  *  ABS  (RN)  -  AINT  (KO  *  ABS  (RN)  /  (2.  *  PI))  * 
k  2.  =*  PI 

TERMl  =  CEXP  (-J  =^  PZQRP)  *  TT 

GO  TO  4 

END  IF 

RMN  =  SQRT  (RN  *  RN  +  ZM  *  ZM  -  2.  *  RN  *  ZM  *  KOS) 

TERMl  =  (0.,0.) 

DO  1  P  --1,  1,2 

DO  1  Q  =  -1,  1.2 

PZQR  =  P  ^  ZM  +  Q  *  RN 
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PZQRP  =  PZQR  *  KO  -  AINT  (PZQR  *  KO  /  (2.  *  PI))  *  2.  *  PI 
TERMl  =  TERMl  +  P  *  Q  *  CEXP  (J  *  PZQRP)  * 
k  EI  (KO  *  (RMN  +  PZQR)) 

1  CONTINUE 

4        TERM2  =  TERM2  +  D  (N)  *  TERMl 

2  CONTINUE 

INTGR  =  INTGR  +  C  (M)  *  TERM2 

3  CONTINUE 

RINTG  =  REAL  (INTGR) 
END 

REAL  FUNCTION  IMINTG  (X,  NX) 
COMPLEX  J,  EI,  TERMl,  TERM2,  INTGR 
INTEGER  NX,  M,  N,  P,  Q 
REAL  *  4  PSIl,  PSI2,  Si,  Tl,  S2,  T2 

REAL  *  4  C  (3),  D  (3),  R,  Z,  RT,  ZT,  U,  V,  CSEC,  COT,  KOS,  Y 
REAL   *  4  X  (NX),  HI,  H2,  Wl,  W2,  KO,  PZQR,  SI,  CI 
&  ,TT,  ZM,  RN,  RMN,  TESC,  TESD,  PI,  PZQRP 

COMMON/  ZINC/  H1,H2,W1,W2,S1,S2,T1,T2.PSI1,PSI2 
COMMON  /PARAM2/  C,  D,  RT,  ZT,  CSEC,  COT,  KOS,  J,  KO 
EI(Y)  =  CI(ABS(Y))-J*SI  (Y) 

u  =  xfn 

V  =  X  (2) 

INTGR  =  (0.,0.) 

TESC  =  ABS  (C  (2)  /  C  (1)) 

TESD  =  ABS  (D(2)  /D(l)) 

TT  =  ALOG  (ABS  ((^l.+KOS)/(l.-KOS))) 

PI  =  4.  *  ATAN  (1.) 

DO  3M  =  1,3 

IF  (M  .EQ.  2  .AND.  TESC  .LE.  l.E-6)  GO  TO  3 

Z  =  (M-1)  *  HI 

ZM  =  -U  *  Wl  *  COT  +  V  *  W2  *  CSEC  +  ZT  +  Z 

TERM2  =  (0.,0.) 

DO  2  N  =  1,  3 

IF  (N  .EQ.  2  .AND.  TESD  .LE.  l.E-6)  GO  TO  2 

R=  (N-1)  *H2 

RN  =  -U  *  Wl  *  CSEC  +  V  *  W2  *  COT  +  RT  +  R 

IF  (ABS  (KO  *  RN)  .LE.  0.5E-1)  THEN 

PZQRP  =  KO  *  ABS  (ZM)  -  AINT  (KO  *  ABS  (ZM)  /  (2.  *  PI))  * 
k  2.  *  PI 

TERMl  =  CEXP  (-J  *  PZQRP)  *  TT 

GO  TO  4 

END  IF 

IF  (ABS  (KO  *  ZM)  .LE.  .5E-1)  THEN 

PZQRP  =  KO  *  ABS  (RN)  -  AINT  (KO  *  ABS  (RN)  /  (2.  *  PI))  * 
k  2.  *  PI 

TERMl  =  CEXP  (-J  *  PZQRP)  *  TT 

GO  TO  4 

END  IF 

RMN  =  SQRT  (RN  *  RN  +  ZM  *  ZM  -  2.  *  RN  *  ZM  *  KOS) 
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TERMl  =  (0.,0.) 
DO  1  P  =  -l,  1,2 
DO  1  Q  =  -l,  1,2 
PZQR  =  P  *  ZM  +  Q  *  RN 

PZQRP  =  KO  *  PZQR  -  AINT  (KO  *  PZQR  /  (2.  *  PI))  *  2.  *  PI 
TERMl  =  TERMl  +  P  *  Q  *  CEXP  (J  *  PZQRP)  * 
k  EI  (KO  *  (RMN  +  PZQR)) 

1  CONTINUE 

4        TERM2  =  TERM2  +  D  (N)  *  TERMl 

2  CONTINUE 

INTGR  =  INTGR  +  C  (M)  *  TERM2 

3  CONTINUE 

IMINTG  =  AIMAG  (INTGR) 

END 
C 

C         Subroutine  to  compute  a  sequence  of  estimates  ESTl  (K)  and 
C  EST2  (K),  1  .LE.  Kl  .LE.  K  .LE.  K2  for  the  N-dimensional  integral 

C  Bl        BN 

C  Int  ...  Int   FUN  (xl,  x2,  ...   xN)  dxl  dx2...  dx3 

C  Al        AN 

C         by  Haber's  method. 

C         Ref:      P.  J.  Davis  and  P.  Rabinowitz,  Methods  of  Numerical 
C  Integratation,  Academic  Press,  1984. 

C 

C         For  each  estimate  ESTl  (K),  two  additional  quantities  ERRl(K) 
C  and  DEVI  (K)  are  computed.  If  the  values  of  DEVI  (K)  do  not 

C  vary  by  more  than  105^  between  consecutive  values  of  K,  then 

C  ERRl  (K)  can  be  taken  as  a  reliable  bound  on  the  difference 

C  between  ESTl  and  the  integral.  A  similar  situation  holds  for 

C  EST2.   DEV2.  and  ERR2   (K).   The  total   number  of  functional 

C  evaluations  is  4  *  (Kl  **  N  +  (Kl  +  1)  **  N  +  ...  +  K2  **  N)  and  K2 

C  should  be  chosen  so  as  to  make  this  may  be  halfed  by  eliminating  the 

C  computation  of  the  EST2  (K).  In  other  situations,  these  values 

C  are  much  better  than  the  ESTl  (K).  A  program  FUNCTION  FUN  (X,   N) 

C         must   be  supplied  bv  the  user  with  X  declared  bv  the  statement 
C  DIMENSION  X  (N).  FUN  must  be  declared  EXTERNAL  in  the  calling 

C  program.   If  N  <  1  or  N  >  10  or  Kl  <  1  or  K2  <  Kl,  the  program 

C         terminates  with  IND  =  0.  Otherwise  IND  =  1. 
C 

C         Modified  by  R.  Janaswamy  so  that  the  output  is  average  of  ESTl 
C  EST2.  Also,  Kl  =  K2  =  K. 

C 

SUBROUTINE  HABER  (N,  LL,  UL,  FUN,  K,  RESULT) 

INTEGER  N,  IND.  KEY.  I,  K,  J 

DOUBLE  PRECISION  AL  (10),  BE  (10),  GA  (10),  B,  G 

REAL  FUN,  Yl,  Y2,  Y3,  Y4,  ESTl,  EST2,  ERRl,  ERR2, 
k  DEVI,  DEV2.  RESULT 

REAL  *  8  SI.  S2,  Dl,  D2 

REAL  LL  (N).  UL  (N),  DEX  (10),  PI  (10),  P2  (10),  P3  (10),  P4  (10), 
k    Ql  (10).  Q2  (10),  Q3  (10),  Q4  (10),  RAN  (10),  AKN,  AK,  T,  JAC 

REAL  AK1,BK 
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EXTERNAL  FUN 

DATA  AL/. 4142135623730950,  .7320508075688773,  .2360679774997897 

.6457513110645906, .3166247903553998,  .6055512754639893. 

.1231056256176605,  .3589989435406736, .7958315233127195, 

.3851648071345040/ 
IND  =  0 

IF  (N  .LT.  1  .OR.  N  .GT.  10)  RETURN 
IND  =  1 
JAC  =  1. 
DO  1  1  =  1,N 
BE  (I)  =  AL  (I) 
GA  (I)  =  AL  (I) 
RAN(I)  =  UL(I)-LL(I) 
JAC  =  JAC  *  RAN  (I) 
DEX  (I)  =  0. 


AK  = 

=  FLOAT  (K) 

KEY 

=  0 

AKl 

=  AK- 

-1.1 

SI  = 

0. 

Dl  = 

0. 

S2  = 

0. 

D2  = 

0. 

AKN 

=  AK 

**N 

T  =  SQRT  (AKN)  *  AK 

BK  =  1.  /  AK 

KEY  =  KEY  +  1 

IF  (KEY  .EQ.  1)  GO  TO  6 

KEY  =  KEY  -  1 

J  =  1 

IF  (DEX  (J)  .GT.  AKl)  GO  TO  8 

DEX  (J)  =  DEX  (J)  +  1. 

GO  TO  6 

DEX  (J)  =  0. 

J  =  J  +  1 

IF  (J  .LE.  N)  GO  TO  4 

GO  TO  3 

DO  71  =  1,N 

B  =  BE(I)  + AL(I) 

IF(B  .GT.  1.)  B  =  B-1. 

G  =  GA  (I)  +  B 

IF(G  .GT.  1.)  G  =  G-1. 

BE  (I)  =  B  +  AL  (I) 

IF  (BE  (I)  .GT.  1.)  BE  (I)  =  BE  (I)  - 

GA  (I)  =  BE  (I)  +  G 

IF  (GA  (I)  .GT.  1.)  GA  (I)  =  GA  (I)  - 

PI  (I)  =  (DEX  (I)  +  G)  *  BK 

Ql  (I)  =  LL  (I)  +  RAN  (I)  *  PI  (I) 

P2  (I)  =  (DEX  (I)  +  l.-G)  *BK 

Q2  (I)  =  LL  (I)  +  RAN  (I)  *  P2  (I) 

P3  (I)  =  (DEX  (I)  +  GA  (I))  *  BK 

Q3  (I)  =  LL  (I)  +  RAN  (I)  *  P3  (I) 
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P4  (I)  =  (DEX  (I)  +  1.  -  GA  (I))  *  BK 
7        Q4  (I)  =  LL  (I)  +  RAN  (I)  *  P4  (I) 
Yl  =  FUN  (Ql.  N) 
Y2  =  FUN  (Q2,  N 
Y3  =  FUN  (Q3,  N 
Y4  =  FUN  (Q4,  N) 

51  =  SI  +  Yl  +  Y2 

Dl  =  Dl  +  (Y1-Y2)**2 

52  =  S2  +  Y3  +  Y4 

D2  =  D2  +  (Yl  +  Y3  -  Y2  -  Y4)  **  2 

GO  TO  5 
3        ESTl  =  0.5  '  SI  /  AKN 

ERRl  =  1.5  ^  DSQRT  (Dl)  /  AKN 

DEVI  =  ERRl  *  T 

EST2  =  0.25  ^  (SI  +  S2)  /  AKN 

ERR2  =  0.75  *  DSQRT  (D2)  /  AKN 
2        DEV2  =  ERR2  -^  T  *  AK 

RESULT  =  0.5  *  (ESTl  +  EST2)  *  ABS  (JAC) 

RETURN 

END 
C 

C  SUBROUTINE  VOLT  WHICH  CALCULATES  THE  VOLTAGE  MATRIX 

C  VM 

C 

C  GEOMETRY  IN  THE  Y-Z  PLANE 

C 
C 

SUBROUTINE  VOLT(F,THITA0,PHI0,A,B,V) 

REALTSI,TCO,PS,PC,KY.KZ,A,B.K0,F,KX 

REAL  MAR.GTUPRS.PRCMOD.KZHTA 

REALS(1:2.30)T(1:230).PSI(1:230),LS(1:230),WI(1:230),THITA0,PHI0 

INTEGER  NMAX 

COMPLEX  V(1:230),J,STR.VOM 

COMMON/  BRAVO/  KX,KY,KZ 

COMMON/  RC2/  LS,WI,S,T,PSI 

COMMON/  RCl/  NMAX 

COSD(X)  =  COS(X='RAD) 

SIND(X)-:SIN(X*RAD) 

PI=4.='ATAN(1.) 

RAD=PI/180. 

K0=PrF/15. 

J=CMPLX(0.J.) 

TSI=SIND(THITAO) 

TCO=COSD(THITAO) 

PS=SIND(PHIO) 

PC=COSD(PHI0) 

KY=KO*TSrPS 

KX=KO*TSrPC 

KZ=K0*TCO 

DO  234  K=1,NMAX 

PRS=SIND(PSI(K)) 
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PRC=COSD(PSI(K)) 

MAR-KY*S(K)+KZ*T(K) 

STR=CEXP(-J*MAR) 

VOM=(STR/SIN(LS(K)*K0))*(A*PRS+B*PRC) 

KZHTA-KY*PRS+KZ*PRC 

MOD=COS(K0*LS(K))-COS(KZHTA*LS(K)) 

GTI=KZHTA**2-K0=^*2 

IF(ABS(KZHTA-K0)/K0.LT.l.E-2.OR. 
k  .ABS(KZHTA+K0)/K0.LT.1.E-2)THEN 

V(K)=VOM*LS(K)*SIN(K0*LS(K)) 

ELSE 

V(K)=(2*K0/GTI)*VOM*MOD 

ENDIF 
234       CONTINUE 
RETURN 
END 
C 

SUBROUTINE  CLUDCP  (A,  N,  NP,  INDX) 

INTEGER  NP,  N,  INDX  (NP),  I,  J,  K,  P,  NMAX 

PARAMETER  (NMAX  =  230) 

COMPLEX  A  (NP,  NP),  TEMP,  ETA,  W  (NMAX) 

REAL  AAMAX,  DUM 

DO  1  K  =  1,  N-1 

AAMAX  =  CABS  (A  (K,  K)) 

P  =  K 

D0  2I  =  K  +  1,N 

DUM  =  CABS  (A  (I,  K)) 

IF  (DUM  .GT.  AAMAX)  THEN 

AAMAX  =  DUM 

P  =  I 

ENDIF 

2  CONTINUE 
INDX  (K)  -  P 
D03  J  =  1,N 
TEMP  =  A  (K,  J) 
A(K,  J)  =  A  (P,  J) 
A  (P,  J)  =  TEMP 

3  CONTINUE 
D04  J  =  K+1,N 
W(J)  =  A(K,  J) 

4  CONTINUE 

D0  5I  =  K+1,N 
ETA  =  A  (I,  K)  /  A  (K,  K) 
A  (I,  K)  =  ETA 
D06  J  =  K+1,N 
A  (I,  J)  =  A  (I,  J) -ETA  *W(J) 
6        CONTINUE 

5  CONTINUE 
1        CONTINUE 

RETURN 
END 
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SUBROUTINE  CLUBSB  (A,  N,  NP,  INDX,  X) 

INTEGER  N,  NP,  I,  J,  INDX  (NP),  NMAX 

PARAMETER  (NMAX  =  230) 

COMPLEX  A  (NP,  NP),  X  (NP),  TEMP,  Y  (NMAX) 
C 

C  DO  PERMUTATIONS  ON  THE  EXCITATION  VECTOR  USING  THE 

C  INFORMATION  ON  THE  ROW  OPERATIONS  DONE  IN  CLUDCP. 

C 

DO  1  K  =  1,  N-1 

TEMP  =  X  (K) 

X(K)  =  X(INDX(K)) 

X  (INDX  (K))  =  TEMP 

1  CONTINUE 
C 

C  FORWARD  ELIMINATION 

C 

DO  21  =  1,N 

Y(I)  =  X(I) 

DO  2  J  =  1,1-1 

Y(I)  =  Y(I)-A(I,J)*Y(J) 

2  CONTINUE 

C 

C  BACK  SUBSTITUTION 

C 

D03I  =  N,  1,-1 

X  (I)  =  Y  (I) 

DO  4  J  =  1+1,  N 

X(I)  =  X  (I) -A  (I,  J)  *X  (J) 
4        CONTINUE 

X(I)  =  X(I)/A(I,I) 

3  CONTINUE 
RETURN 
END 

C 

C  SUBROUTINE  RADIAT  TO  COMPUTE  THE  RADAR  CROSS  SECTION 

C  OF  A  PLANAR  STRUCTURE  COMPOSED  OF  ARBITRARILY 

C         ORIENTED  PLANAR  DIPOLES. 

C 

SUBROUTINE  RADIAT  (F,THITAO,PHIO,NMAX,A,B,RCS) 

REALRC,THITA0.PHI0.LS(1:230).WI(1:230),S(1:230),T(1:230),P1 

REALCR,F,K0.PI.UN,DM,EM,RAD,P,R,A,B,THITA,PHI,RCS,PSI(1:230) 

REAL  KX.KY,KZ.KRON,KG 

INTEGER  G.NMAX 

COMPLEX  J,NTHI0,NPHI0,TEMP1,TEMP2,CUR(1:230) 

COMMON/  BRAVO/  KX,KY,KZ 

COMMON/  RC2/  LS,WI,S,T,PSI 

COMMON/  RC4/  CUR 

COSD(X)  =  COS(X*RAD) 

SIND(X)  =  SIN(X*RAD) 
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PI  =  4.  *  ATAN(1.) 

RAD  =  PI/180. 

KO  =  PrF/15. 

UN=120.*PI 

THITA=180.-THITA0 

PHI  =  180.  +  PHIO 

J=CMPLX(0.,1.) 

TEMP1  =  (0.,0.) 

TEMP2=(0.,0.) 

NTHIO=(0.,0.) 

NPHI0=(0.,0.) 

D0  58G=1,NMAX 

P1=PSI(G) 

DM=K0*(S(G)*SIND(THITA)*SIND(PHI)+T(G)*COSD(THITA)) 

EM=K0*(SIND(P1)*SIND(THITA)*SIND(PHI)+ 
&  +C0SD(P1)^SIND(THITA)) 

P=SIND(P1)*C0SD(THITA)*SIND(PHI)-C0SD(P1)*SIND(THITA) 

TEMP1-P*CEXP(J*DM)*CUR(G)/SIN(K0*LS(G)) 

TEMP2=CEXP(J*DM)*CUR(G)*SIND(P1)*COSD(PHI)/SIN(KO*LS(G)) 

IF(ABS(EM-K0)/K0.LT.l.E-2.OR.ABS(EM+K0)/K0.LT.l.E-2)THEN 

KG-LS(G)*SIN(KO*LS(G)) 

ELSE 

KG=-2*K0*(COS(EM*LS(G))-COS(K0*LS(G)))/(EM**2-K0*"2) 

ENDIF 

NTHI0=TEMP1*KG+NTHI0 

NPHI0-NPHI0+KG*TEMP2 
58       CONTINUE 

CR=CABS(NTHI0)**2+CABS(NPHI0)**2 
KRON=(A^KY+B*KZ)/KX 

RC=(K0**2)*(UN**2)*CR/(4.*Pr(ABS(A)**2+ABS(B)**2+ 
k  +ABS(KR0N)**2)) 

IF((RC/((30./F)**2)).LT.l.E-6)THEN 

RCS  =  10.*ALOGlO(l.E-6) 
ELSE 

RCS=10.*ALOG10(RC/((30./F)**2)) 
ENDIF 
RETURN 
END 
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APPENDIX  B 
PROGRAM  TREE 


C  PROGRAM  TREE 

C 

C  PROGRAM  TO  GENERATE  A  PLANAR  FRACTAL  TREE  WITH 

C  REDUCTION  FACTOR  (COND)  AND  BRANCHING  ANGLE  THETA 

C  WRITTEN  BY  T.R.  NELSON,  PhD,  UNIVERSITY  OF  CALIFORNIA 

C  SAN  DIEGO,  LA  JOLLA,  CA,  90293,  AND  MODIFIED  BY 

C  LT.  JOHN  DEMIRIS  TO  GENERATE  THE  INPUT  DATA  FOR 

C  THE  RCS  PROGRAM. 

C 

C  INPUT  DATA: 

C 

C  L  =  NUMBER  OF  BRANCHING  LEVELS. 

C  N  =  NUMBER  OF  BRANCH  SEGMENTS. 

C  ILEN  =  INITIAL  LENGTH. 

C  ISP  =  INITIAL  STARTING  POINT. 

C  COND  =  REDUCTION  FACTOR. 

C  THETA  =  BRANCHING  ANGLE. 

C 

C  OUTPUT  DATA: 

C 

C  IX,  lY  =  COORDINATES  OF  FIRST  AND  END  POINT  OF  EACH 

C  GENERATED  BRANCH. 

C  INPUT  DATA  FOR  RCS  PROGRAM 


C 


REAL  ZZ  (1:5,1:1024),ISP,ILEN,IX,IY,ITH1 

OPEN  (UNIT=1.  FILE  =  'INGEOM',FORM  =  'FORMATTED') 

OPEN  (UNIT=2,  FILE  =  'OUT',  FORM  =  'FORMATTED') 

OPEN  (UNIT=10,FILE  ='INPUT1',F0RM='F0RMATTED') 

READ(1,^END=10)L,N,ILEN,ISP,COND,THETA 

CLOSE(l) 

ZZ(  1,1)  =  ISP 

ZZ  (4,1    =  ILEN 

ZZ  (3,1)  =  0. 

ZL  =  ZZ(4,1) 

ZZ  (2,1)  =  0. 

ZN  =  N 

PI  =  3.14159265 

IX  =  0. 

lY  =  O.-ILEN 

IX  =  ZZ(2,1) 

IY  =  ZZ(1,1) 

DO  100  LOOP  =1,L 
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ZLOOP  =  LOOP 

NPTS  =  N**ZLOOP 

INC^O 

NPREV  =  N**(ZL00P-1) 

PREV  =  NPREV 

ANGLE  =  (ZN-L0)/2.0 

ZINDX  =  -ANGLE 

DO  200  J  =1,NPTS 

ZJ- J 

IF  (INC.GE.N)  ZINDX=-ANGLE 

IF  (INC.GE.N)  INC=0 

IW=NPTS-J+1 

IR  =  PREV-ZJ/ZN+1.0 

ZL=ZZ(4,IR) 

T  =  ZZ(3,IR) 

X  =  ZZ  1,IR 
Y  =  ZZ(2,IR) 
ZLl  =  ZL*COND 
IF(ZLLLT.O.l)  GOTO  300 
THl  =  ZLl/33.0 

Tl  =  T+THETA*ZINDX 

ZZ(3,IW)  =  Tl 

T2  =  T1*PI/180.0 

IX  =  Y 

IY  =  X 

TEMPX1=IX 

TEMPY1=IY 
C 

C  IX,IY  ARE  THE  COORDINATES  OF  THE  FIRST  POINT  OF  THE 

C  BRANCH  ALONG  THE  X,Y  AXIS  RESPECTIVELY 

C 

WRITE(2,*)  IX,IY 

XI  =  ZL1*C0S(T2)+X 
Yl  =  ZL1*SIN(T2)+Y 
IX=Y1 

IY=X1 
C 

C  IX,IY  ARE  THE  COORDINATES  OF  THE  END  POINT  OF  THE 

C  BRANCH  ALONG  THE  X,Y  AXIS  RESPECTIVELY 

C 

WRITE(2,*)  IX,IY 

TEMPX2=:IX 

TEMPY2=IY 
C 

C  S  AND  T  ARE  THE  COORDINATES  OF  THE  CENTER  OF  EACH 

C  BRANCH 

C 

S  =  (TEMPX2+TEMPXl)/2. 

T  =  (TEMPY2+TEMPYl)/2. 

PSI1=ATAN((TEMPX2-TEMPX1)/(TEMPY2-TEMPY1)) 

PSI=PSI1*180./PI 


ZLHALF=ZLl/2. 

THlHALF=THl/2 

WRITEllO.^')  ZLHALF,TH1HALF,S,T,PSI 

ZZ(1.I\V)=X1 

ZZ  2J\V)=Y1 

ZZ(4J\V)-ZL1 

ZINDX=ZINDX+1.0 

INC=INC+1 
200  CONTINUE 

100      CONTINUE 
300      STOP 
END 
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